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(57) Abstract 



A method is described for estimating the bone quality of a vertebrate, on the basis of two-dimensional image data comprising 
information relating to the trabecular structure of at least a pan of a bone of the vertebrate, the image data being data obtained, by exposing 
at least the pan of the bone to electromagnetic radiation, the method comprising subjecting the image data to a statistical analysis. 
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A METHOD OF ESTIMATION 

The present invention relates to the field of non- invasive 
determination of bone quality of vertebrates. This determina- 
tion may be used in the diagnosis of, e.g., osteoporosis and 
5 other bone diseases which cause bone fragility and, thus, 
increase the risk of bone fracture. 

Accelerated bone loss leading to osteoporosis is a common 
phenomenon in women after menopause. Women often ignorantly 
suffer from accelerated bone loss, and the reduced strength 

10 of the bones is not discovered until a bone is broken or a 
vertebra collapses due to a load which a healthy bone or 
vertebra should be ahle to withstand. Thus, a large part of 
the osteoporotic patients are not aware of the reduced 
strength of their bones until a fracture reveals the disease 

15 and the extent thereof . 

Today, osteoporosis affects more than one third of elderly 
women in the industrialized part of the world. The prevalence 
of this disease is still increasing, partly caused by the 
increase in the proportion of elderly people, partly for 

20 unexplained reasons. In the lesser developed parts of the 
world, like South-east Asia and South America, it is pre- 
dicted that osteoporosis will become an enormous socio-econo- 
mic burden within the next 20-30 years. However, if indivi- 
duals at risk of developing osteoporosis can be identified, 

25 preventive measures can be applied. This requires a reliable, 
cheap and safe method for identification of those at risk. 

Since the risk of osteoporosis is closely related to the bone 
resistance to fracture, bone strength or bone quality mea- 
surements are likely to supply the essential information. 

30 Since the means of preventing osteoporosis are much more 

efficient than those of treating osteoporosis, identification 
of individuals at risk of developing this disease is crucial. 
Only by early prevention of osteoporosis, the individual as 
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well as the socio-economic consequences of the disease may be 
minimized . 

The method according to the invention provides a method for 
detection of reductions in bone quality, such as architecture 
5 and strength, at an early stage. Thus, the method of the 

invention provides a tool of screening potential patients for 
bone diseases and, thus, provides the possibility of early 
detection of bone diseases at a stage in which no symptoms of 
disease are noted by the patient . 

10 A typical measure of bone strength has been the Bone Mineral 
Density (BMD) of the bone. HMD measurements are typically 
obtained by X-ray of a bone together with a standard wedge. 
Having determined which part of the wedge attenuates the X- 
ray beam to the same degree as the illuminated bone, a 

15 measure of BMD may be obtained. 

A more sophisticated method of determining BMD is by X-ray 
absorptiometry of two different wavelengths . Using two wave- 
lengths encOales the method to compensate for the effects of 
soft tissue etc. aroxind the bone and, thus, to obtain a more 
exact determination of the X-ray attenuation of the bone. 



20 



However, as will be clear from the following, the Bone Mine- 
ral Density of a bone is not necessarily connected to the 
actual strength of the bone. The reason for this is to be 
found in the structure of a bone of a vertebrate. 

25 A bone of a vertebrate consists of a cortical outer layer and 
a ccmcellous inner structure. Omission of the cancellous 
inner structure of a bone would result in a quite fragile 
bone. The cancellous inner structure of the bone consists of 
so-called trsibeculae. Thicker vertical trabeculae are 

30 positioned in the bone in the direction of the main load 
(main compression or pull in the bone) , and thinner, 
horizontal trabeculae interconnect the vertical trabeculae. 
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Thus, the main density of a bone is constituted by the corti- 
cal layer and the vertical trabeculae. This is quite natural, 
as the bone reserves the largest part of the material to 
withstand the most common loads . A BMD measurement 
5 determining the density of a bone will therefore primarily 
estimate the amount of bone in the cortical layer and the 
vertical trabeculae and, thus, only to some degree the 
ability of the bone to withstand the loads to which the bone 
is adapted. 

10 However, when loads are applied not in the direction of the 
vertical trabeculae, the structure of the bone would be 
fragile without the thinner, horizontal trabeculae which 
interconnect the vertical traQjeculae. The horizontal 
trabeculae define the ability of the bone to withstand loads 

15 not in the direction of the vertical trabeculae . Thus , the 
strength of the bone in directions not in the direction of 
the vertical trabeculae is defined by a very small part of 
the bone mass. 

Different bones vary in the relation between cortical and 
20 trabecular bone structures. Furthermore,, the architecture of 
the trabeculae vary according to the potential loading of 
this part of the bone. The number and diameter of 
longitudinal trabeculae auid the presence of transverse 
trcQseculae are also of fundamental importance for the 
25 fragility of the bone. These differences are of major 

importance for the strength of the bone, but they can only to 
a minor extent be detected by BMD measurements . 

Since trabecular bone has a greater surface area/weight 
ratio, it is to a higher extend exposed to accelerated bone 

30 loss than the cortical bone. Therefore, osteoporotic frac- 
tures are mainly seen in bones with predominantly trabecular 
stiructure , such as the distal forearm, proximal humerus, 
thoracic and lumbar spine as well as the femoral neck. 
Changes in BMD related to loss of trabecular bone severely 

35 underestimates the loss of resistance to fracture. As a 
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consequence, BMD measurements of patients with osteoporotic 
fractures reveal a major overlap with BMD measurements of 
normal subjects of the same cohort. 

The cibove information explains why persons with low energy 
5 fractures may have normal BMD values when con5)ared to age- 
matched normal controls. Thus, one may wish to distinguish 
between bone quantity and bone quality, where the bone qual- 
ity is more closely related to the mechanical or bio-mechan- 
ical strength of the bone . 

10 A measure of the overall bone strength may naturally be 
obtained from a bone specimen taken from the potential 
patient and subjected to mechanical testing. However, this 
rec[uires bone biopsy, which is painful and implies a minor 
risk of complications. Thus, in order to have a comfortable, 

15 cheap, fast and safe screening of the very large group of 

potential patients (most women after menopause) , the estima- 
tion of the bone quality should be performed non-invasively 
in a safe and fast manner. 

The present invention offers a non- invasive method which 
20 provides measures of the overall strength of bones and which 
is safe, comfortcdsle for the potential patient, fast and 
cheap . 

In "A New Method for Automatic Recognition of the Radio- 
graphic Trabecular Pattern" Wil G.M. Geraets et al . , Journal 

25 of Bone and Mineral Research, Vol. 5, No. 3, 1990, pp 227- 
233, a method of recognizing the trabecular pattern from X- 
ray pictures is disclosed. According to this method, two 
types of noise are removed from an X-ray picture which has 
been scanned into an image memory: high-frequency noise is 

30 removed by using a median filter, and low-frequency noise is 
removed by using a local averaging operation on the image. 

After noise reduction a segmentation of the image data is 
performed whereby the image is binarized and subsequently 
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eroded (by using a median axis transform) so as to retain 
only central lines having a thickness of 1 pixel . A total of 
7 features are derived from the segmented and the eroded 
images . These features are subsequently correlated to the 
5 Bone Mineral Density measurements of the bones . 

However, due to the almost insignificant correlation between 
the image features obtained and the BMD measurements in the 
above-mentioned reference, it is evident that a higher cor- 
relation and, thus, a more precise estimate of the bone 
10 quality is required. 

In "cotrputerized Radiographic Analysis of Osteoporosis: 
Preliminary Evaluation", Plilip Caligiuri et al , Radiology 
1993, 186 ,-417-474, a method in which X-ray images of the 
lumbar spine are scanned analyzed. In this reference, the 
power spectrum was obtained and the Root Mean Square (RMS) 
and the First Moment (FMO) thereof were con^ared to typical 
BMD measurements performed on the bones , It was concluded 
that both BMD did not correlate well with FMO and RMS and 
that this may be due to FMO and RMS containing additional 
information. In this reference, only the power spectrum and 
features derivable therefrom were investigated. 

The present invention provides a method and a system for 
providing an easily performed, but highly reliable estimation 
of the bone quality. 

In a first aspect, the present invention relates to a method 
for estimating the bone quality of a vertebrate, on the basis 
of two-dimensional image data comprising information relating 
to the trabecular structure of at least a part of a bone of 
the vertebrate, the image data being data obtained by expos- 
ing at least the part of the bone to electromagnetic radi- 
ation, the method comprising subjecting the image data to a 
statistical analysis comprising 
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a background correction procedure in which low frequent 
intensity variations not related to the trabecular struc- 
ture of the bone is reduced relcitive to image data 
related to the trabecular structure of the part of the 
5 bone , 

an image manipulation and feature extraction procedure 
wherein at least the local image intensity information as 
well as variation in the local intensity are utilized to 
extract information related to the trabecular structure 
10 of the part of the bone, the image manipulation and 

feature extraction procedure comprising siibjecting the 
image data to at least one of the following procedures : 

(a) obtaining an estimate of the parametric estimate 
of the power spectarum of the image data and 

15 extracting features relating to the energy 

distribution of the parametric estimate, 

(b) obtaining, on the basis of image data on which a 
Fourier method has been used to en^hasize the 
information in the image data relating to the 

20 trabecular structure, an estimate of a grey- level 

co-occurrence matrix and extracting at least one 
feature on the basis of the estimated co-occur- 
rence matrix, 

(c) obtaining an estimate of the projected trabecular 
25 pattern of the image data by using a Fourier 

method to emphasize the information in the image 
data relating to the trabecular structure and 
subjecting the manipulated image data to a mor- 
phological operation, and extracting features 
30 relating to the trabecular structure from the 

estimated projected trcUaecular pattern, 

(d) obtaining, on the basis of a frequency analysis 
of the image data, features relating to the 
periodicity of the trabecular structure of the 

35 part of the bone, 

and an estimation procedure in which the bone quality of 
the vertebrate is estimated on the basis of the derived 
features and optionally other features related to the 
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bone or the vertebrate and a predetermined relationship 
between the values of such features and reference bone 
quality parameters. 

As indicated above, in the present context "bone quality" is 
5 not equalled to "bone quantity" , such as Bone Mineral Den- 
sity, as even a small loss of bone quantity may lead to a 
large loss of bone quality if the bone loss has taken place 
at critical positions in the trabecular structure. Thus, in 
the present context "bone quality" is a measure closely 
10 related to the risk of fracture of the bone, as it has been 
demonstrated in stress/strain evaluations of bone biopsies in 
which the microscopic structure is also evaluated (See, e.g., 
Lis Mosekilde relating to this issue) . 

Even though the bones and, thus, the trabecular structures 
15 are inherently three-dimensional, a projection of this struc- 
ture into two dimensions, such as a radiographic image, 
yielding a so-called projected trsUaecular pattern, conveys 
sufficient information about the trabecular structure and, 
thus, the bone quality to give a satisfactory estimate of the 
20 bone quality. 

Information relating to the trabecular structure may be local 
variations in, e.g., grey level information or any other 
information from which features relating to the trabecular 
structure may be derived. 

25 Low frequency intensity variations in the image data will 

naturally depend on how the image data are obtained. If the 
image data are obtained on the basis of radiographic images, 
the low frequency intensity variations may be due to, e.g., 
scattering of the electromagnetic radiation, anatomic struc- 

3 0 tures surrounding the illuminated part of the trabeculae such 
as cortical bone, fat tissue and muscles of varying thickness 
as well as, e.g., inhomogeneous X-ray illumination of the 
part of the bone. 
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According to the first aspect of the present invention, local 
image intensity information and variation in the local inten- 
sity are utilized to extract information relating to the 
trabecular structure of the part of the bone. In, e.g., 
5 digitized radiographic images, local image intensity informa- 
tion may be the individual pixel values, whereas variation in 
the local intensity is related to the textural" information 
contained in, e.g., inhomogeneities in the image data. 

The extracted features resulting from the image manipulation 
10 and feature extraction procedure quantify properties of the 
trabecular structure and, thus, of the bone quality. The 
extracted features are subsequently introduced into an esti- 
mation procedure in which a predetermined relationship 
between features and bone quality enables .the estimation 
15 procedure to estimate the bone quality of the vertebrate. 

In the present context, to "emphasize" information, such as 
magnitude information, means to give prominence to prevailing 
frequency information. This may be performed by either enhan- 
cing the prominent information or by reducing the less domi- 
20 nant information - optionally both. 

According to the first aspect of the invention, an estimate 
of the bone quality is obtained in an estimation procedure on 
the basis of a predetermined relationship between features 
obtained and reference bone quality parameters. This prede- 
25 termined relationship is typically estsQslished through stat- 
istical modelling, where explanatory variables (image fea- 
tures and optionally other explanatory features relating to 
the bone quality) are used to model corresponding reference 
bone quality data (response variable) . 

30 As described eOsove, even though bone quality is not equal to 
bone quantity in the present context, the reference bone 
quality parameters are preferably parameters related to the 
strength of the bone rather than, e.g., BMD or BMC which, as 
explained above, relate to the density of the bone rather 



) 
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than the strength thereof. However, to provide even higher 
diagnostic sensitivity, combination of the present invention 
with a measurement of BMC/BMD is contemplated to be useful . 
BMC/BMD measurements can be obtained either by the above - 
5 mentioned two-wavelength- technique or they may, if the image 
data are obtained on the basis of an X-ray image, be avail- 
able from the same image data, where, e.g., a standard alu- 
minum wedge has been illuminated together with the bone in 
question . 

10 At present, the most common way of obtaining non- invasive 

information about the bone quantity and bone quality is based 
on X-ray illumination. X-ray images may easily be used to 
generate the image data for use in the present invention. As 
it is preferred to perform the statistical analysis on a 

15 computer, the analogue X-ray image may be digitized and 
introduced into a computer by scanning the X-ray image. 

Naturally, the resolution of the X-ray film used to obtain 
the X-ray image will have an effect on the quality of the 
image data. Thus, it is presently preferred that the resol- 
20 ution of the X-ray film is at least 4 pairs of lines per 

centimetre, such as at least 5 pairs of lines per centimetre. 
Even though this may be sufficient for obtaining a satisfac- 
tory estimation of the bone quality, it is preferred that the 
resolution of the X-ray film is at least 10 pairs of lines 

25 per centimetre, such as at least 25 pairs of lines per centi- 
metre, prefereibly at least 50 pairs of lines per centimetre. 
It is possible to obtain resolutions of the X-ray film up to 
at least 100 pairs of lines per centimetre, such as at least 
250 pairs of lines per centimetre and probedaly as high as 500 

3 0 pairs of lines per centimetre, such as at least 600 pairs of 
lines per centimetre. Naturally, a resolution of this 
preferred size will ensure that a large amount of information 
is present in the X-ray image. 

Naturally, it is preferred that the scanning of X-ray images 
35 is performed with a sufficiently large resolution in order to 
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ensure that a minimum relevant information of the X-ray image 
is lost in the transfer to the image data. Thus, it is pre- 
ferred that the scanning has been performed at a resolution 
of at least 10 lines per cm, such as at least 25 lines per 
5 cm, preferably at least 100 lines per cm, such as at least 
200 lines per cm, such as at least 250 lines per cm. 

Furthermore, it is preferred that the resolution of the 
scanner is better than 4 true bits per pixel, such as better 
than 6 true bits per pixel, preferably better than 8 bits per 
10 pixel. 

Naturally, in order to have the scanner actually scan the 
image, it should be able to trans illuminate the radiographs. 

The ESKOSCAN 2450 from Eskofot A/S has been found to fulfil 
the above criterias and to be highly suited for use in the 
15 method of the invention. 



It is presently preferred that the background correction 
procedure comprises at least reducing or optionally removing 
low frequency information having a frequency significantly 
lower than the spacing of the projected trabeculae. As 

20 described above, this low frequent spectral content of the 
image data is typically caused by cortical bone, fat tissue 
and muscles of varying thickness as well as, e.g., 
inhomogeneous X-ray illumination of the part of the bone. 
Naturally, this type of undesirsUale effect is unavoidable 

25 when obtaining non- invasive image data on the basis of X-ray 
radiographs. However, other non- invasive image acquisition 
techniques, such as, e.g., MR cind CT imaging may generate 
image data not automatically or not to the same degree "suf- 
fering from" this type of undesired effect. 

30 Even though it is difficult to exactly quantify the limit 

between undesired low frequency information and the desired 
higher frequency relevant information independently of the 
image resolution, it is presently preferred that information 
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having frequencies half or less than the spacing of the 
projected trabeculae is at least reduced or optionally 
removed. More preferably, information having frequencies 
being a quarter or less, such as a tenth or less than the 
5 spacing of the projected trabeculae is at least reduced or 
optionally removed. 

One preferred way of reducing or preferably removing the low 
frequency information is using a background correction pro- 
cedure comprising generating secondary image data as a result 
10 of performing a median filtering with a predetermined kernel ' 
size and subtracting this result from the original image 
data. One of the advantages of using a median filter is that 
this operation is edge preseirving. 

Another way of reducing or preferably removing the low fre- 
quency information is using a background correction procedure 
comprising generating secondary image data as a result of 
performing a mean filtering with a predetermined kernel size 
and subtracting this result from the original image data. A 
mean filtering is typically much faster than the median 
filtering. However, the mean filtering is not edge preser- 
ving. Edges generated in this part of the process may have an 
adverse effect and generate false or erroneous information 
later in the process, leading to wrong estimates of the bone 
quality. 

25 It is typically preferred that the kernel size is at the most 
1/2 of the image data, such as at the most 1/4 of the image 
data, preferably at the most l/lO of the image data, more 
preferadDly at the most 1/20 of the image data. 

A third way of reducing or preferaOaly removing the low fre- 
30 quency information is using a background correction procedure 
comprising globally fitting a two-dimensional polynomial to 
the image data and generating background corrected image data 
on the basis of the residuals of the fitting procedure. Apart 
from potential difficulties in determining the optimal order 



15 
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of the polynomial, which the person skilled in the art will 
know, this method potentially offers an extremely fast back- 
ground correction of image data. 

It is contemplated that a suitable order of the polynomial 
5 may be at the most 15, such as at the most 10, more preferab- 
ly at the most 5 . 

Having performed a background correction of the image data, 
an image manipulation and feature extraction procedure is 
performed in order to obtain features quantifying textural 
10 properties of the trabecular structure. 

An important textural property relating to the bone quality 
is the degree of anisotropy of the projected trabecular 
pattern. This property may be described in an intuitively 
appealing way using the power spectrum of the image data. The 
15 power spectrum may be obtained using direct methods, auto- 
covariam.ce methods or parametric methods . 

Classical spectral estimation, however, which comprises 
direct methods and auto- covariance methods, typically gives 
relatively non-consistent spectral estimates with a poor 
20 resolution. Alternatively, smoothing methods may be employed. 
This, however, brings about the well-known 
variance/resolution trade-off. 

Therefore, parametric spectral estimation is used according 
to the invention as this estimator is a consistent spectral 
25 estimator with superior resolution properties. 

One way of obtaining the parametric estimate of the power 
spectrum is subjecting the image data, optionally weighted 
with a window, to a Fast Fourier Transformation. This is a 
so-called direct method of which there are several variants, 
30 where the most popular probably is the so-called Welch 
method. 
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Another way of obtaining the parametric estimate of the power 
spectrum is subjecting an estimate of the auto-covariance 
function of the image data, optionally weighted with a 
window, to a Past Fourier Transformation. 

5 As mentioned above, the parametric methods potentially offer 
a higher resolution than the classical methods. The basic 
principle of parametric methods is to identify an appropriate 
model of the image data, which is assumed to be a homogeneous 
texture, and subsequently estimate the parameters of the 
10 model. Having estimated these parameters, the parametric 
estimate of the power spectrum of the fitted model may be 
obtained. In one dimension, this type of spectral estimate 
has proven itself superior to the classical spectral 
estimates . 

15 The methods, on the basis of which the parametric estimate is 
obtained are at present preferably chosen from the group 
consisting of: causal Simultaneous Auto -Regressive Moving 
Average (SARMA) models, non-causal SARMA models, Gaussian 
Markov Random Field models and Maximum Entropy Spectral 

20 Estimates. 

When considering spatial data, such as image data, causality 
is not natural. Abandoning causality, however, leads to 
serious problems in the estimation of the parameters as the 
Least Squares estimator is no longer consistent. Therefore, 
25 it is often preferred to superitt^jose an artificial 

directionality on the image data. Thus, having a causal 
model, the Least Square estimation principle may be applied. 

However, a disadvantage of superimposing an artificial 
directionality on the image data is the so-called "direc- 
30 tional bias" which may be quite conspicuous in, e.g., the 
estimate of the parametric estimate of the power spectrum. 
Consequently, the use of causal models in the present context 
is not an obvious choice. This problem is partly solved using 
the "parallel-resistor" averaged spectral estimator 
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introduced by Chan in ("Two-dimensional Spectral estimation 
from Auto-regressive Models with varying Areas of Support" 
from 1981) . 

Hitherto, when considering non-causal models, one has been 
5 forced to consider the Maximum Likelihood Estimator, which 
gives a number of practical problems for the following rea- 
sons : 



test for this in two dimensions or more, due to the 
computational effort required. 

For the above-mentioned reasons, the Maximum Likelihood 
estimation is not a practical way of obtaining estimates of 
25 the model. Other estimation principles may, of course, be 
more practical alternatives . 

At present, it is preferred that the method on the basis of 
which the parametric spectrum is obtained, is a non- causal 
SARMA model using the so-called MORSE estimator. This estima- 
30 tor will be described further below. 



10 



15 



20 



The Likelihood function includes the determinant of the 
Jacobian matrix which is of the size M^xM^ (where M in 
the 2-D situation is the side length of the image) . This 
determinant is a complicated non-linear function of the 
model parameters . 

Even if assumptions are made about the image (e.g. torus) 
that offers an analytical expression for the above-men- 
tioned determinant, the computational load will still be 
very large. 

Close to non-stationarity, the Likelihood function is 
extremely rippled. This will often cause the non- linear 
optimization routine to find a sub-optimal solution to 
the estimation problem. It is not feasible/possible to 



Having obtained the parametric estimate of the power spectrum 
of the image data, one or more features are extracted which 
quantify properties of the bone quality. 
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The distribution of the energy in the parametric estimate of 
the power spectrum relates directly to the anisotropy of the 
image data. Thus, several types of features may be proposed 
which quantify this anisotropy. 

5 At present, it is preferred _to obtain at least one feature 
related to parameters of a contour encompassing at least a 
predetermined percentage of the energy of the parametric 
estimate of the power spectrum. 

This the contour may, e.g., is determined by 

10 (a) defining a contour, in terms of one axis in each 

dimension, around the centre of the parametric 
estimate of the power spectrum, all points on the 
contour having substantially the same distance to the 
centre of the parametric estimate of the power 

15 spectrum, and the contour encompassing less than the 

predetermined percentage of the energy of the 
parametric estimate of the power spectrum, 

(b) for each dimension of the data calculating the per- 
centage of the energy of the parametric estimate of 

20 the power spectrum encompassed by a dilated contour 

in which the axis in the dimension in question is 
increased by a predetermined distance, 

(c) increasing the axis of the contour, in the dimension 
in which the increase in energy encompassed by the 

25 dilated contour is the largest, by the predetermined 

distance and 

(d) repeating steps (b) and (c) until the percentage 
encompassed, by the contour exceeds or equals the 
predetermined percentage. 



30 



The distribution of energy in the parametric estimate of the 
power spectrum is typically oblong. Tests have shown that the 
distribution of the energy in the parametric estimate of the 
power spectra of osteoporotic patients differ significantly 
from those of non-osteoporotic persons. This difference may. 
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e.g., be seen in a difference in the area covered by or the 
elongatedness of the aibove- mentioned contour. 



One way of obtaining suitsUble features from the contour is to 
define the shape of the contour and subsequently derive 
5 features from parameters of the contour. 

A suitaJDle way of defining the contour is to define the axes 
of the contour so that these are orthogonal and directed 
along principal directions of the parametric estimate of the 
power spectrum. 

10 A natural shape of this contour may be an ellipsoidal contour 
and one or more features may be derived from the semi -axes 
thereof. However, a rectangular contour may also be used. 

The above-mentioned method of defining the contour is prefer- 
ahly performed so as to obtain the smallest possible contour 
15 of the desired shape. This will ensure that the contour is 
uniquely defined and that the features derived from parame- 
ters of the contour relate as well as possible to the anisot- 
ropy of the image data. 



An alternative way of extracting information in the frequency 
20 domain is to derive information from e.g. the height, width 
or total area under peaks of a smoothed periodogram. From 
this information, features relating to periodicity of the 
trabeculae may be derived. In this manner, information 
relating to a periodicity of the trabeculae is obtained by 
25 performing a frequency analysis of the image data. 



Several other ways of extracting features from the parametric 
estimate of the power spectrum may be used, such as fitting a 
Gaussian to the normalized parametric estimate of the power 
spectrum and deriving features from this Gaussian. Other 
30 methods are using rings and wedges as described by Weszka et 
al. (1976) . 
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An image manipulation and feature extraction method alterna- 
tive to or supplementing the above use of the parametric 
estimate of the power spectrum is to estimate the projected 
trabecular pattern of the image data by using a Fourier 
5 method to emphasize the information in the image data 
- relating to the trabecular structure "and siibj ecting the 
manipulated image data to a morphological operation, and 
extracting features relating to the trabecular structure from 
the estimated projected trabecular pattern. 

10 The emphasizing procedure is essential in order to bring out 
the projected trabecular pattern in a robust manner. To this 
end, Fourier methods have proven suitable, and it is pre- 
ferred that the Fourier method comprises 

subjecting the image data to Fourier transformation, 
15 - performing a subsequent mathematical transformation of 
the Fourier transformed data, 

converting the transformed data back into the spatial 
domain . 

Naturally, the subsequent mathematical transformation may be 
20 linear or non-linear. 

Even though it is preferred that the subsequent mathematical 
transformation is performed by raising the magnitude of the 
Fourier transformed data to a power larger than 1, such as a 
power in the range of 1.1-10, preferably 1.5-4, such as sibout 
25 2.0, also other transformations, such as other types of 
filtering, may be performed. 

In order to further emphasize the information relating to the 
trabecular staructure, a third mathematical transformation may 
be performed in order to preserve only part of the magnitude 
30 information of the transformed data. This may, e.g., be 

removal of at least part of the information in the magnitude 
of the spectrum. 
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Having emphasized the inforroation in the image data relating 
to the trabecular structure, the resulting image data is 
preferably subjected to grey- scale morphological operations 
in order to prepare the data for feature extraction. 

5 At present, the purpose of the grey- scale morphological 
operations is to extract a binary representation of the 
projected trabecular pattern (PTP) . The foreground colour 
represents the trabeculae, whereas the background colour 
represents the cavities. 

10 The above PTP is preferably obtained using the so-called 

morphological top-hat operation followed by a thresholding. 
The result of these operations is preferably further 
"cleaned" by removing small isolated clusters of pixels. 

On the basis of the PTP obtained, a number of different 
15 features may be obtained. In the following, the generation of 
just a few of the vast number of possible features is 
described. 

One method of generating features from the PTP is performing 
an operation (the so-called "distance transform") comprising 
20 determining, for each background pixel, the distance accord- 
ing to a given metric from the background pixel to the 
nearest foreground pixel. The features extracted from this 
operation are indeed intuitively appealing as they relate 
directly to the inter- trabecular distance. 

25 The above metric is presently the Euclidian distance. How- 
ever, other metrics may eventually prove more suitable in the 
present context . 

Features based on the disteuice transform may, e.g., be 
derived from a mean value and/or standard deviation and/or 
30 the coefficient of variation and/or skewness and/or kurtosis 
of the determined distances. 
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Another method of generating features from the PTP is per- 
forming an operation (hereinafter referred to as "star area 
transform" ) comprising determining a measure for each back- 
ground pixel based upon determining the distance from the 
5 background pixel to the nearest foreground pixel in a number 
of given directions in the image data.' 

Again, features derivable from the star area transform may, 
e.g., relate to a mean value and/or standard deviation and/or 
the coefficient of variation and/or skevmess and/or kurtosis, 
10 optionally the maximum distance, of the determined measures. 

In addition to the above-mentioned features derivable from 
the image data, it may be preferred to input additional data 
into the estimation procedure. It is evident that a large 
number of features which relate to the illuminated bone, but 
15 which are not derivable from the image data, may enhance the 
precision of the estimation of the bone quality. 

As the present method may be used to estimate the bone qual- 
ity of any bone from any vertebrate, additional information 
as to the age and/or sex, and/or species, and/or race and/or 
20 the specific bone considered in the vertebrate, and/or a 
estimated Bone Mineral Density, and/or a estimated Bone 
Mineral Content, may be included in the estimation procedure. 

Even though the BMD may be introduced in the estimation 
procedure, this measure may also optionally be determined 
25 from the image data . BMD may be estimated by including data 
from a reference object in the exposure of the bone to the 
electromagnetic radiation and on the basis of the absorption 
of the electromagnetic radiation of the bone and of the 
reference object. 

30 The purpose of the estimation procedure is to estimate given 
bio-mechcuiical properties of the bone on the basis of the 
introduced features extracted from the image data and other 
explanatory variaJsles . 
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One way of obtaining an estimation procedure of the above 
type is to have the estimation procedure based on a statisti- 
cal model , taking into account the correlation structure in 
the data set, so as to assign appropriate weights to the 
5 significant features in accordance with the predetermined 
relationship . 

The cLbove- mentioned model may be determined in a number of 
ways. However, estadDlishing a model of this type requires 
obtaining corresponding values of all relevant features and a 
10 response variable relating to the bio-mechanical property 
which is to be estimated by the estimation procedure. The 
bio-mechanical property in question may be an absolute or a 
relative measure of the bone quality. 

Examples of interesting bio-mechanical properties of bone are 
the mechanical bone strength, and/or a Bone Mineral Density 
measurement, and/or a Bone Mineral Content measurement, 
and/or a score value by a skilled radiologist. Naturally, the 
method of the present invention will output a value corre- 
sponding to the calibration of the method: if the method is 
calibrated towards a strength parameter, the output of the 
method will relate to the strength of the bone in question. 

An important purpose of the present invention is to provide a 
method for evaluating bone strength and fracture risk of 
significant clinical value, thus to provide a higher estima- 
25 tion of fracture risk than the best technique availc±)le 

today, the Double X-ray Absorptiometry Bone Mineral Content 
measurements (DXA-BMC measurements) . 

As the calibration of the estimation procedure will determine 
the estimated parameter of the bone, this parameter should be 
30 chosen in a mcumer so that it correlates to the bone quality. 
The assessment of the extent to which the embodiments of the 
present invention are eible to estimate bone quality may be 
performed by: 
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1. Cohort- studies , in which occurrence of new fractures are 
recorded and related to the initial estimation. 



2. cross -sectional studies, in which fracture cases and 
age- and- sex-matched controls are compared to estimations 

5 obtained by the method according to the invention and option- 
ally further to DXA-BMC measurements, or 

3. by using bio-mechanical testing of bone specimens, pre- 
ceded by estimating the bone quality according to the inven- 
tion and optionally additionally DXA-BMC-measurements of the 

10 bone in question. 



Using either distal radius, lumbar spine, femoral neck, or 
any bone in which two-dimensional image data comprising 
information relating to the trabecular structure the bone can 
be obtained, coherent values for determined bone strength, 
15 DXA-BMC and the estimate according to the present invention 
may be performed on the basis of human post-mortem bone 
samples or samples from other vertebrates . 

The calibration of the method of the present invention may 
also be calibrated along the above- illustrated methods as 
20 these methods generate bone quality information and estima- 
tion on the basis of corresponding image data. 

Bone strength may, e.g., be evaluated by directly measuring 
the plasticity and maucimum load of a bone in a stress -strain 
diagram (See, e.g.. Lis Mosekilde) . These parameters are 
25 preferably obtained in both the direction of the vertical 
trabeculae and in the orthogonal direction in order to have 
the most complete estimate of the bone quality. 

In the method of the invention, sxibstantially all types of 
models may be used in the estimation procedure. At present, a 
30 preferred model is selected from the group consisting of: a 
General Linear Model , a Generalized Linear Model, an Artifi- 
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cial Neural Network, a Causal Probabilistic Net or Classifi- 
cation And Regression Trees. 

As described above, the present method may be used to esti- 
mate the bone quality of bone in any vertebrate. Thus, the 
5 vertebrate may be a human, a horse, a great ape, a large ape, 
an anthropoid ape, a pig, a cow, etc, and the actual bone may 
naturally be virtually any bone in the vertebrate, such as 
radius, femur, corpus vertebrae, calcaneus, talus, os carpi, 
metatarsi, metacarpi, falanges, tibia, fibula, patella, ulna, 
10 humerus, mandible, clavicula, scapula, os coxae, os 

naviculare, os cuboideum, os cuneiform I, os cuneiform II or 
OS cuneiform III. 

A second aspect of the present invention relates to a method 
for estimating the bone quality of a vertebrate, on the basis 

15 of two-dimensional image data comprising information relating 
to the trabecular structure of at least a part of a bone of 
the vertebrate, the image data being data obtained by expos- 
ing at least the part of the bone to electromagnetic radi- 
ation, the method comprising subjecting the image data to a 

20 statistical analysis comprising 



25 



a background correction procedure in which low frequent 
intensity variations not related to the trabecular struc- 
ture of the bone is reduced relative to image data 
related to the trabecular structure of the part of the 
bone, 

an image manipulation and feature extraction procedure 
comprising subjecting the image data to at least one of 
the following procedures: 

(a) obtaining an estimate of the parametric estimate 



30 



of the power spectrum of the image data and 
extracting features relating to the energy 
distribution of the parametric estimate, 



35 



(b) obtaining, on the basis of image data on which a 

Fourier method has been used to emphasize the 
information in the image data relating to the 
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trsUDecular structure, an estimate of a grey- level 
co-occurrence matrix and extracting at least one 
feature on the basis of the estimated co-occur- 
rence matrix, 

5 (c) obtaining an estimate of the projected tralaecular 

pattern of the image data by using a Fourier 
method to emphasize the information in the image 
data relating to the trabecular structure and 
subjecting the manipulated image data to a mor- 
10 phological operation, and extracting features 

relating to the trabecular structure from the 
estimated projected trabecular pattern, 
(d) obtaining, on the basis of a frequency analysis 

of the image data, features relating to the 
15 periodicity of the trabecular structure of the 

part of the bone, 
and an estimation procedure in which the bone quality of 
the vertebrate is estimated with a Multiple Correlation 
Coefficient better than 0.5 on the basis of the derived 
20 features and optionally other features related to the 

bone or the vertebrate and a predetermined relationship 
between the features and reference bone quality parame- 
ters . 



In fact , the bone quality of the vertebrate is preferably 
25 estimated with a Multiple Correlation Coefficient better than 
0.55, such as better than 0.6. Naturally, the higher the 
Multiple Correlation Coefficient the higher the correlation 
between the actual bone quality and the estimated bone qual- 
ity. Thus, a Multiple Correlation Coefficient better than 
30 0.65, such as better than 0.7, or better than 0.8, such as 
better than 0.85 is preferred in order to have as high a 
correlation to the actual value as possible. 

In a third aspect, the present invention relates to a method 
for estimating Non-Causal Simultaneous Auto- Regressive Moving 
35 Average models in two or more dimensions, the method compris- 
ing optimizing a given direct measure of the flatness of the 
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residual spectrum of the model. In the present context 
"optimizing a given direct measure of the flatness" means 
that the flatness is maximized. 

Non-linear optimization procedures are typically used to 
maximize measures of the above type. A suitable non-linear 
procedure may comprise the following steps 

(a) generating a set of initial parameters for the model, 

(b) generating the residual spectrum of the model on the 
basis of the parameters, 

(c) obtaining the measure of the flatness of the residual 
spectrum, 

(d) obtaining a new iterate of the parameters on the 
basis of the flatness measure and a search direction 
in parameter space, 

(e) repeating steps (b) - (d) until given stop criterion is 
reached . 

A number of measures of the flatness of the residual 
parametric estimate of the power spectrum may be chosen. A 
well known measure of flatness of probaJaility distributions 
20 is the so-called entropy. Therefore, in the present context, 
the entropy is measured on the normalized residual parametric 
estimate of the power spectrum, disregarding the DC value. 

Preferred embodiments of aspects of the invention will now be 
described with reference to the drawings wherein 

25 Fig. 1 shows a typical radiographic image of a wrist, 
obtained by an ESKOFOT scanner (See Example l) , 

Fig. 2 shows an extracted sub- image of the image of Fig. l, 

Fig. 3 shows the image of Fig. 2 after a median filtering 
with a 25x25 kernel size. 



10 



30 Fig. 4 shows the final background corrected image, 
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Fig, 5 shows the power spectrum of the background corrected 
image of Fig. 4, 

Fig. 6 shows the restored result of the image of Fig. 4 after 
having been sxibjected to an emphasizing procedure, 

5 Fig. 7 shows the result of a morphological top-hat operation 
on the restored image of Fig. 6, 

Fig. 8 shows the image of Fig. 1, wherein small isolated 
groups of pixels (noise) have been removed (this image is 
referred to as the projected trabecular pattern (PTP) ) , 

10 Fig. 9 shows the result of a distance transform calculated 
for the PTP of Fig. 8, 

Fig. 10 illustrates the star area transform for a single 
pixel . 

Fig. 11 shows the result of the star area transform of all 
15 background pixels in the PTP of Fig. 8, 

Fig. 12 shows the result of the maximum distance transform of 
the PTP of Fig. 8, 

Fig. 13 shows the smallest possible ellipse containing 70% of 
the spectral energy for a 'normal' patient, 

20 Fig. 14 shows the smallest possible ellipse containing 70% of 
the spectral energy for an osteoporotic patient. 

Fig. 15 illustrates the systematic bias of the LS estimates 
for simulated isotropic textures (the vertical lines along 
the curve outline the empirical confidence intervals for the 
25 mean of the estimated parameters { , 

Fig. 16 illustrates the systematic bias of the approximative 
ML estimates (close to non- stationarity) for simulated 
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isotropic textures (the vertical lines along the curve out- 
line the empirical confidence intervals for the mean of the 
estimated parameters) , 

Fig. 17 illustrates a part of the likelihood function for a 
5 first-order bilateral SAR model (given a realization of the 
process with the parameters a=/3=-0 . 22735 . The likelihood is 
plotted in the interval (a, jS) € [ -0 . 275 , -0 . 20] x [-0 . 275 , -0 . 20] ) , 

Fig. 18 illustrates that good estimates may be obtained using 
the torus-ML estimator with the approximate ML estimates as 
10 initial values combined with a stationarity check (if a 
non- stationary solution model is found, the optimizer is 
restarted with new initial values; the vertical lines along 
the curve outline the empirical confidence interval for the 
mean of the estimated parameters) , 

15 Fig. 19 illustrates that good estimates may be obtained using 
the MORSE estimator (the vertical lines along the curve is 
the empirical 95% confidence intervals for the mean of the 
estimated parameters) , 

Fig. 20 shows the objective functional of the MORSE estimator 
20 close to non- stationarity (as opposed to, e.g. , the ML-estim- 
ator, this objective functional is not rippled, which makes 
it much easier to perform the non-linear optimization) , 

Fig. 21 illustrates the periodogram estimator of the power 
spectrum for a normal individual (left) and for an 
25 osteoporotic individual (right) . The spectrum is shown in 

inverse so that the high intensity values are darker than the 
lower intensity values, 

Fig. 22 illustrates the parametric NSHP estimate of the power 
spectrum for a normal individual (left) and for an 
30 osteoporotic individual (right) . The spectrum is shown in 

inverse so that the high intensity values are darker than the 
lower intensity values, 
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Fig. 23 illustrates the parametric WQP estimate of the power 
spectrum for a normal individual (left) and for an 
osteoporotic individual (right) . The spectrum is shown in 
inverse so that the high intensity values are darker than the 
5 lower intensity values, 

Fig. 24 summarizes an experiment performed in order to 
describe the applicability of the MORSE estimator, 

Fig. 25 illustrates histograms based on parameters estimated 
for each of 1000 synthesized 128x128 textures, 

10 Fig. 26 illustrates optimal strengths plotted against the 
predicted strength, where age is not included in the model, 

Fig. 27 illustrates optimal strengths plotted against the 
predicted strength, where age is included in the model, and 

Fig. 28 illustrates the correlation of the logarithmic 
15 optimal fracture load and a single feature suitably chosen. 

EXAMPLE 1: Textural analysis o£ bone structiure 

In the following Example, a preferred embodiment of the 
method of the invention will be described wherein X-ray 
images of the distal radius and a vertebral body from the 

20 lumbar spine (L3) are restored using image processing. The 
restored images are used as a basis for extracting textural 
features that are shown to correlate well with the structure, 
the density and the fracture load of the trabecular bone. 
Several ways of extracting promising textural features are 

25 described. A model based on textural features (and possibly 
age and sex) is described and shown to discriminate well 
between osteoporotic and non-osteoporotic cases. 
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Obtaining suitable imaae data 

The X-ray images used in this Example were obtained using a 
stsindard X-ray instrument having a focus area of 0.6 ram^, a 
tube voltage of 45 kV, a focus-to-foil distance of 100 cm, a 
5 single, fine X-ray foil (it should be noted that no foil is 
required) , using a double emulsion film (it should be noted 
that a single emulsion film is equally applicsUole) and adjus- 
ting the MAS product to obtain a suitably exposed film. 

The X-ray pictures were scanned into a computer using a 
10 ESKOFOT ESKOSCAN 2450 scanner capable of scanning up to 5000 
lines/cm. However, at present 600 lines/cm was used. This 
scanner was found highly suited for the present use. A typi- 
cal image is shown in Fig. 1. 



The demands to a scanner suitable for this use are that it 
15 should naturally be able to transilluminate the X-ray film, 

it should preferably have a resolution better than the grains 
in the film/foil combination (so as to loose no information 
in the scanning procedure) and it should preferably generate 
grey- values of at least 8 true bits, in order to obtain a 
20 suitable dynamic range of the image data. 



Spatial Texture Models 



Preprocessing 



In this study, only the textural information is taken into 
account. Therefore a sub- image is extracted in a region of 

25 the bone containing only tredsecular bone, as e.g. shown in 
Fig. 2. Low frequent intensity variations due to, e.g., 
scattering of the radiation, anatomic structures as cortical 
bone, muscles, fat tissue of varying thickness etc., are 
removed by subtracting a median filtered version of the image 

30 (shown in Fig. 3) from the image itself (a technique known as 
"unsharp filtering") . In this case a 31x31 median filter is 
used. The background- corrected image is shown in Fig. 4. This 
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image is very noisy and should be further processed prior to 
further transformations and subsequent feature extraction. 

The noise in the image is clearly very high-f requent , whereas 
the information which is preferably considered (the trabecul- 
5 " ar structure) is of a more' low-frequent nature. 

The power spectrum of the extracted, background corrected 
sub-image is shown in Fig. 5. In the following, an emphasiz- 
ing procedure will be described, which may not be generally 
applicable for removal of high-frequent noise, but tends to 
10 emphasize the dominating (in this case low- frequent) features 
in the images, while suppressing less dominating (in this 
case high frequent noise) features. This property makes it 
attractive for this specific purpose. 

Algorithm 1 The following steps outline the restoration 

15 algorithm that removes high-f requent noise and emphasize the 
relevant structures in the image . 

1. Perform a 2-D FFT and obtain the phase, Xp(fj., f c) < 
the magnitude, Xf^itj., f^) • 

2 . The magnitude is raised to the power p : 



20 3. The histogram of S(fr, fc) is stretched to match a Gaus- 
sian: SQ(fj-, fc) 

4. A context image, C{t^, t^) . is formed by thresholding 
SQ(fj., fc> preserving T% of the pixels. Objects less than n 
pixels are removed from the thresholded image, (note that 

25 more than T% of the energy is preserved) . 

5. A 'new' magnitude image. XN(fr. ^c) ' is formed using the 
context image and the power spectrum (not the magnitude) : 
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S{f^.f^) if C{f^.f^)=l 



0.0 otherwise 



(2) 



6. The new magnitude image X^Cf^, f^) is combined with the 
phase image Xp(fj., f^), and the restored image, XR(r,c), is 
obtained as the real part of the inverse Fourier transform. 
It is preferred to use thresholding in the range T=l-6 and 
5 p=l-2. 

Note that the magnitude, raised to the power p, is used to 
form the new magnitude image. This may not be obvious, but it 
has been experienced that this in general provides very- 
useful results. The reason for this is, that when the 
10 magnitude is raised to a given power larger than one, the 

parts of the spectrum that contain most of the energy tend to 
be emphasized whereas spectral components with less energy 
tend to be weakened further. The result of the enphasizing 
procedure using T=2 and p=2 is shown in Fig. 6. 

15 

Extraction of the PTP 

In the following, various image transformations and textural 
features, extracted from transformed images, will be dis- 
cussed. Several features are based upon an image subsequently 
2 0 referred to as the projected trcOaecular pattern (PTP) which 
is obtained as outlined below: 

Algorithm 2 The steps outlined below describe the desired 
transformations of the restored image in order to measure a 
relevant set of features on the resulting PTP. 

25 1. Perform the grey-scale morphological top-hat opera- 



tion on the restored image. The top-hat operation is 
the original image minus an opened version of the 
image. This image is subsequently thresholded. The 
result of these operations is shewn in Fig. 7. 
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2. This image clearly contains a lot of noise-pixels. 
Small isolated objects (in the present Example all 
objects less than 10 pixels) are removed. The result- 
ing PTP is shovm in Fig . 8 . 

5 Extraction of Statistical Features 

Naturally, it is preferred to extract features that, in some 
sense, quantify the properties (such as structure and 
density) of the trabecular structure. Not surprisingly, the 
meshes of the PTP-net are much larger, and often broken, for • 
10 osteopenic individuals than for non-osteoporotic individuals. 
This difference may be quantified in several ways. Such 
properties are, e.g., reflected in the average distance from 
each background pixel to the nearest foreground pixel 
(bone-pixel) in the projected treibecular pattern. 

15 The size of the dark meshes in the pattern is a relevant 

parameter that is conveniently quantified by calculating the 
so-called distance transform for the projected trabecular 
pattern. The result of the distance transform on the PTP is 
shown in Fig. 9 and is denoted xd^^'*^) ■ 

20 In addition, for each background pixel in the projected 
treibecular pattern, a line is drawn from the pixel, in a 
given direction, until a foreground pixel (' bone -pixel ' ) is 
found. The distance between these pixels is measured. This is 
done for a fixed number of directions as illustrated in Fig. 

25 10. After this, the average distance which is referred to as 
the star area, denoted Xp^i^.c) , or the maxintum star length 
which is denoted Xx^^'^^ ^® calculated. These images are 

shown in Fig. 11 and Fig. 12 respectively. 

A number of statistical moments may be calculated from the 
30 resulting images from the measured average, standard 

deviation, coefficient of variation and skewness for all 
background pixels: 
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(3) 



z, c 



d. = 



(4) 



a. 



cC. 



(5) 



r (x.(r,c) 



(6) 



where the siibscript denotes different transforms considered: 

• D: Distance transform 

• M: Mean value of star length 

• S : Standard deviation of star length 
5 • X: Maximum star length 

• .A: Star area 

We briefly outline the interpretation of the considered 
parameters : 

The [i. and o^. features are either related to the general mesh 
10 size (for D, M, A and X) or the degradation of the mesh (S) . 
For example, larger mesh sizes are reflected in a larger 
value of ft.. Thus, for increasing values of JT. suid o^. , the 
texture index should decrease . 

The cv. feature is somewhat different. It measures the width 
15 of the eirpirical distribution in the transformed images. If 
all the distances are large, cv. would be small. The same 
thing would be the case, if all the distances were small. 
Typically, however, we find that the empirical distribution 
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is wide, only if the FTP mesh is broken. Thus, for increasing 
values of cv. , the texture index should typically decrease. 

The 7. feature is also related to the mesh size. If the 
histogram of distances is heavier in the right side, 7.<0. If 
- 5 the empirical distribution is symmetrical, 9=0- For large 
meshes in the PTP, the distribution of distances would 
typically be skewed for large values. Thus, for increasing 
values of j". , the texture index should decrease . 

Tffytrural Features Exfrracted in the Spatial Frequency Domain 

10 One would intuitively think that features measured in the 

frequency domain should carry relevant information about the 
trabecular structure. The trabecular structure for strongly 
osteoporotic patients is clearly more low frequent than that 
for 'normal' patients. In the present section, the considered 

15 Fourier features are not measured on the projected trabecular 
pattern. Instead, the features are measured on the spectrum 
of the background corrected image. 

Preprocessing 

Before estimating the spectrum, non-stationaries of the 
20 extracted image are removed by subtracting a 31x31 median 

filtered version of the image from itself as described aibove. 

Spectral Estimation 

As the main interest is directed toward the relatively low- 
frequent properties of the power spectrum, a weighed estimate 
25 of the dispersion may be used, where p*(u,v) is the 

normalized power spectrum taking a weight function into 
account : 
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p'iu, V) = 



P(u, v) • w{u, v) 



<U,v)«(0,0) > ' ' \ t I 



(7) 



and 



ur( u, v) = — -i 



(8) 



Given a positive definite symmetrical matrix 



E- 



/ 2 \ 



(9) 



the eigenvalues are: 
^~ 2 



(10) 



_ (ou-^gt-v/(o5-ot)'-^4oQ 

' 2 



(11) 



5 and the corresponding (normalized) eigenvectors 



p7^ 



(12) 



An alternative estimate of the power spectrum may be obtained 
from the following process: 



The probably best known estimator of the power spectirum is 
the periodogram estimator: 



wo 96/07161 PCT/DK95/00338 

35 



J(u,v) = ^ 



M-l M-1 



m-0 n-o \ / 



(13) 



where (u,v) are spatial frequencies. 

It is well known that this estimate is a non- consistent 
estimate of the power spectrum. For example, when the samples 
from a Gaussian white noise process are used, the periodogram 
5 estimate has a variance that is proportional to the variance 
of the white noise process . 

One approach to variance reduction is to segment the data, 
such as in non- overlapping blocks, where the spectrum is 
estimated in each block, using this periodogram estimator, 
10 and obtaining the final estimate by averaging the 

periodograms of the blocks. In this manner, the variance is 
decreased by a factor equal to the number of blocks . Of 
course, this is at the expense of frequency resolution. 

In the following, an iterative process is suggested for 
15 estimating the dispersion matrix, which estimate is more 

sensitive to the anisotropies reflected in the low frequency 
parts of the orientation of the texture. 



Algorithm 3 The dispersion matrix of the approximating 

Gaussian is estimated iteratively in the following way: 

20 1. = I 

2 . i .= 1 

3 . repeat 

3a. Estimate the dispersion: 



I 2 ^ 



2 



(14) 



where 



) 
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U V 
U V 

^uv=5252 '^^ (17) 

U V 

and 



(u. V) =(u,v)g^"j (18) 



3b. i = i-1 
until I |Ei - Ei.J < 6 

From both these techniques, the principal direction of the 
5 power spectrum may be found using: 



a = cos 



-1 ( Ouv ] 



This principal direction may be used to find the ellipse with 
the smallest possible area, containing a fixed fraction of 
the energy of the original (non- weighted) power spectrum. The 
ellipse, rotated with the angle a around the origin, is 
10 expressed in polar coordinates: 



i?(9)2 = (20) 

iJ^cos^ (e-o) +a2sin2 (e-a) 



where a and b are the semi-major and semi -minor axis, 
respectively, and where the above-mentioned X's are the 
eigenvectors . 



In the present Example, the semi -major and semi -minor axes 
15 may be found using the following scheme: 



I I 
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Algorltimi 4 The steps outlined below are used to estimate 

the semi-major and semi-minor axis in the ellipse with the 
smallest possible area, containing a given fraction of the 
spectral energy. The angle with which this ellipse should be 
5 rotated around the origin is found as outlined above. 

1. Start with the semi -minor and semi -major axis suffi- 
ciently small, e.g., a = b = 10 . 



2. Increase a and b with one, and calculate the fraction 
of energy in each ellipse (rotated with the angle a 

10 around the origin) . 

3 . The ellipse with the largest fraction of energy is 
chosen. 



4. Repeat the steps 2 and 3 until the requested fraction 
of energy is reached. 

15 

The features considered here are the estimated values: a, b, 
a-b, b/a as well as the features X^, X2, X2/X1, P=o^^^/a^a^, 

cos-l( ^""^ 1 (21) 



-if . °- . ] 



(o^-ot)'-4ot. 



(22) 



535^ (u2 + vr2)p(u, V) (23) 

U V 

-5252p(u, V) log(p(u, v) ) (24) 

U V 



In Fig. 13, the estimated ellipse containing 70% of the 
20 energy (auid the direction) is shown for a 'normal' patient. 
In Fig. 14, the 70% ellipse is shown for an osteoporotic 
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patient. Clearly, the frequency content for the osteoporotic 
patient is a lot more low-frequent. This means that the 
energy is concentrated near the centre of the image and, 
consequently, the areas of the ellipses are very different 
although they contain the same fraction of energy. 

Spatial Interaction Mndf^l^ 

Above, it is described that a powerful set of features may be 
extracted in the spatial frequency domain. However, the 
performance of features extracted in the frequency domain 
depend heavily on the spectral estimator used. The 
periodogram is a non-consistent spectral estimator. This 
problem may be overcome by segmenting the data in non- 
overlapping segments, estimating the periodogram in each 
segment and finally averaging over these periodograms . 

15 The variance reduction is directly proportional to the number 
of segments. However, this is at the expense of frequency 
resolution . 

In the following, we will consider two types of causal 
Simultaneous AutoRegressive (SAR) models to obtain parametric 
20 spectral estimates,- the Non- Symmetrical Halfplane (NSHP) 
model and the Quarter-Plane (QP) model. 

Assume a wide sense stationary, zero mean, stochastic process 
defined on a rectangular grid of pixel sites D = {s= (s^., s^.) -. 
OsSj-sM, Oss^sN} , obeying a SARMA model defined as: 



3'-* E *ry,-r= E (25) 

25 where r={rj.,r^) e N^r or and Uf/^j^ are the regions of 

support for the parameter arrays {4>} and [6], respectively. 
We also refer to these sets a neighborsets which, in some 
sense, defines the order of the model. 
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The innovation process £g is assumed to be an IID N(0,Eg^) or 



Gaussian white noise sequence. If = 0, we have a pure SAR 

model, and if N;^ =0, we have a pure SMA model. 

In the following, we will focus on pure autoregressions , for 
5 the following reasons: The use of mixed SARMA models would 
probably give more parsimonious models, but the models 
required to capture the nature of the considered texture are 
of a such low order that parsimony in not really an important 
issue. It also turns out that the important feature presently 
10 is the shape of the spectrum in terms of e.g. peak locations." 
Also, the LS estimator for the pure SAR model is a set of 
linear equations, whereas non- linear optimization is required 
in order to obtain estimates of mixed SARMA models . 

The (non -normalized) spectral density function is given as 
15 the squared modulus of the frequency response function 
multiplied by the noise variance: 




(26) 



2 {i^ei£:iN^) ,rHQ^2iN^) r)^ 



where the spatial frequencies (fj.,fj,) = feCif, and 



l__r 

2 M' 2 M 



) : iz, c) eQ 



(27) 



The vector notation is defined as: 



(28) 



(29) 
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£(N.) f=£:£2l[cO3(0i[f.r) ) ,r€N.] 



(30) 



£i(N. ) f=£:al[sin{(ji {f, r) ) ,t€N.] 



(31) 



« (f , r) =2it {r^'ff*r^-f^) 



(32) 



So far, the question of causality has not been addressed. In 
two dimensions, the notion of causality is not natural. 
Abandoning causality, however, leads to a number of 
difficulties, one of them being that the Least Squares 
5 estimator of the parameters is no longer consistent. Hence, 
it has become popular to impose an artificial directionality 
on the image, thus obtaining a causal model (in fact, 
recursively computatcdale is a more correct term) . 

The two types of causal models typically considered in the 
10 literature are the Quarter Plane support models (QP) and the 
Non-Symmetrical Half Plane support models (NSHP) . In fig. 21, 
parametric estimates using the NSHP model are illustrated. It 
is seen that the spectrum is slightly elongated in a certain 
direction, when compared to the corresponding periodogram 
15 estimate seen in Fig. 22. This phenomenon is due to 

directional bias which is an effect that is even more 
pronounced for the QP mode. In order to solve this problem, a 
Weighted QP (WQP) estimator has been proposed: 



where the subscripts ++ and -+ simply refer to the considered 
20 quadrant, i.e. the first and the second quadrant. This 

parametric estimator is widely considered to be among the 
best parametrical spectral estimators, but it suffers from 
the drawback that it has no model interpretation. 

In fig. 23, a WPQ estimate is illustrated in which no 
25 directional bias is seen. 



2 




(33) 
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Using these estitnates, features of the type described above, 

that is: X^, Xj, X2/X1, p=aj.^/a^a^, a^/Oj?, and: 



o=cos'^ 



(34) 



(35) 



(fr-f'o)p(fr.fc) (36) 

-EE^<^r'-fc)log(p(f,,f,)) (37) 
may be used. 

Also, as OL describes the principal orientation of the 
5 approximating Gaussian, this direction may be used to find 
the cLbove -mentioned energy ellipses. Thus, also the semi- 
major and semi-minor atxes of these ellipses S and 6, 
respectively as well as the features a-b, b/a may be used. 

Extraction of features relating to periodicity 

10 An alternative method of deriving features from the 

preprocessed image of Pig. 6 is to again subject this image 
to the two dimensional Fourier transformation and extract the 
Fourier power spectrum (the periodogram spectrum) . 

It is presently preferred that the resolution of the image is 
15 reduced using a Gaussian pyramid, thus, preserving as much of 
the information as possible in the high resolution prior to 
the Fourier transform. 

Feature Eictractioji 



It is not surprising that the power spectrum comprises peaks 
20 corresponding to periodicity of the vertical trabeculae. It 
is, however, fortunate that also pea)cs - however small - 
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corresponding to a periodicity of the horizontal trabeculae 
may be discerned from the power spectrum. 

Information relating to the density of trabeculae, their 
thickness and the distance between them may be extracted from 
5 the height of these peaks relating to the periodicities, the 
total area under these peaks and/or the "width" thereof . 

At present, it is preferred not to extract the peak positions 
from the Fourier power spectrum, as the periodogram or the 
power spectrum are very noisy due to the fact that the 
10 periodogram estimator of the power spectrum is non- 
consistent. Instead, a combined algorithm is preferred, where 
the parametric spectrum is used to identify these peak 
positions, and the smoothed periodogram is used to estimate 
the associated energy etc. 

15 In the power spectrum, a total of four peaks are typically 
present - or more precisely, only two peaks, as the 
parametric spectrum is symmetric. A line is drawn through the 
two largest peaks (which relate to the vertical trabeculae) 
through the DC-value. From the middle of. this line, a line is 

20 drawn perpendicularly thereto also through the DC-value. It 
is inherent from the model that the peaks of the horizontal 
trabeculae are positioned perpendicularly to those of the 
vertical traibeculae. The msucima on the new line will 
correspond to the peaks of the horizontal trabeculae. 

25 This method is preferred in order to be independent and 
insensitive to any rotation of the image. 

Smoothing of the spectrum may be performed using a 3cxk 
kernel, such as a 5x5 kernel using, e.g., a mean filter, a 
median filter, a Gaussian filter etc. 

30 From this spectrum, the peak heights, the volume under the 
peaks a.s.o. may be determined and used in the prediction of 
the bone quality. 
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In addition, it is contemplated that also the curvature of 
the peaks in the parametric power spectrum may be used as a 
predictor of bone quality, as this curvature also relates to 
deviations of the periodicity. 

5 EZAMPIiE 2 : Study of. feature perf ormazice 

A data set was used to investigate different image features 
and compare their performance . Radiographic images were 
obtained from 97 individuals and each image was ranked by 
experienced radiologists. The analysis of this data set is 
10 described in the following. 

Each patient was assigned the score +1 if the radiologist 
considered the patient to be 'normal' and the score -i in 
case of accelerated bone loss (osteoporotic) . The four radio- 
logists agreed upon the diagnosis in 52 cases, and only these 

15 are used in this preliminary study. In this study of feature 
performance, a general linear model with the score as 
response variable and age, sex and textural features as 
explanatory variables is used. Furthermore, rigorous statis- 
tical tests of the significance of the individual features is 

2 0 performed using a backwards elimination procedure. 

The features used have been described above, and are all 
based on first, second and third order statistics measured on 
transforms of the projected trabecular pattern (PTP) or 
fitted ellipses on the Fourier power spectrum of the restored 
25 images . 

The fitted models are of the type 

Si = Po * Pi^g^e^ + Pz^e-^i * P3^i,i "* * €i (3 8) 

where s^ is the score of the i'th patient, age^ the age of 
the i'th patient, sex^ the sex of the i'th patient (+1: Male, 
-1: Female in order to make both sexes contribute equally to 
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the mean) and fj,i the j ' th image feature of the i'th 
patient. It is assumed that c is i.i.d. gaussian white noise. 
Under this assumption, the LS estimates are equal to the ML 
estimates. Test statistics are calculated for each parameter 
5 in the model. Using a backwards elimination approach, the 
least significant parameter of the model is eliminated, and 
the model is subsequently re-estimated with the remaining 
parameters. This is continued until no more parameters can be 
eliminated from the model. 

10 Features based on the projected trabecu lar nattif^.rn 

Out of 12 image features, 9 were eliminated using the back- 
wards elimination procedure. The order in which they were 
eliminated is summarized in the table below: 



Feature 
number 


Feature 
name 


1 


2 




4 


o 




1 


8 


9 


I" 


1 








1 










2 


od 








1 • 












CVd 








1 












4 


Id 






• 














5 


Us 




t 


1 


• 








6* 
























CVs 


• 








1 








8 


7s 




1 


• 




1 








9 


Hm 












• 






10 












1 






• 


11 






• 






1 1. 








12 


-f.\t 


i" ■ 






1 






• 





15 The boldfaced features marked with an asterisk are signifi- 
cant parameters that could not be eliminated. 

The resulting model is given by 



A summary of the resulting model is given below: 
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Call: liii( formula = score [idx4] - 

age + sex[idx4] + f[idx4, l] + f [1(3x4, 3] 
+ f [iclx4, 6] ) 



Residuals 



Min IQ Median 3Q Max. 

-3.755 -0.8549 -0.06046 0.7199 3.854 



Coefficients : 





Value 


Std. Error 


t value 


Pr(>| t| ) 


( Intercept) 


50 .2650 


12 .2399 


4 .1067 


0.0002 


age 


-0 .0995 


0 .0141 


-7 .0427 


0 .0000 


sex [idx4] 


1 . 1937 


0.2579 


4 . 6281 


0 . 0000 


f[idx4, 1] 


-11.7739 


2 .3699 


-4.9681 


0 .0000 


f[idx4, 3] 


-37 .2419 


11.0728 


-3 .3634 


0 .0016 


f[idx4, 6] 


2 .7775 


1.2370 


2 .2454 


0 .0297 



15 Residual standard error: 1.746 on 45 degrees of freedom 
Multiple R-Squared: 0.7516 

Correlation of Coefficients: 

(Intercept) age sex(idx4] f[idx4,l] f[idx4,31 

age 0.0063 

20 sex[idx4] -0.1329 0.2430 

f[idx4, 1] -0.6991 -0.1573 -0.0231 

f[idx4, 3] -0.9837 -0.0538 0.1220 0.7305 

ftidx4, 61 -0.0287 0.1205 0.1072 -0.6269 -0.0918 



Most of the explanatory effect is contained in the age and 
25 the sex. The rest of the variability (approximately 16% 
point) is explained by the textural features. It may be 
desired to transform some of the features, or even derive new 
ones . 



From the cibove, it may be concluded that the considered 
30 statistics measured on at least the distance transform and 
the star area of the projected trabecular pattern, are cjuite 
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Strong features. This investigation, however, does not sug- 
gest the maximum distance transform as an optimal feature. 



Fourier Features 



In the following, four features extracted from the power 
5 spectrum will be considered. The only feature eliminated from 
this set of features is the length of the semi-major axis 
(a) . This is in good agreement with what would be expected, 
as the erosion of the trabecular structure in osteoporotic 
patients is mainly in the vertical direction. The performance 
10 of the resulting model is comparable to the results obtained 
using features based on the distance transform. 

Again the data were analyzed in the statistical software 
package Splus, and the following output was obtained for the 
resulting model : 



15 Residuals : 





Min 


IQ Median 3Q 


Max 






-2 . 852 


-1.114 -0. 


3153 0.6817 


4.854 






Coefficients 














Value 


Std. Error 


t value 


Pr (>| t 1 ) 


20 


(Intercept) 


8 .1977 


2 .9152 


2 .8120 


0 .0073 




ages [idx4 , 2] 


-0.1086 


0.0153 


-7 .1036 


0.0000 




sex [idx4] 


1.2674 


0 .2711 


4 .6747 


0.0000 




f req tidx4 , 2] 


8 . 6433 


2 .5788 


3 .3516 


0.0016 




f req [idx4 ,3] 


-0 .2102 


0 .0660 


-3 .1845 


0.0026 


25 


f req [idx4 ,4] 


-88 . 6948 


24.0649 


-3 .6857 


0.0006 



Residual standard error: 1.843 on 45 degrees of freedom 
Multiple R-Squared: 0.7232 
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Correlation of Coefficients: 

(InurcEpt) •gc*|i<lx4.2| 
0 .0165 

0.0103 0.2150 

-0.2471 -0.1272 

0.2189 0.1303 

0.1145 0.086b 



ages [ idx4 , 2] 
sex [ idx4 ] 
freq[idx4,2l 
freqlidx4, 31 
freq(idx4, 4] 



tu(lilx4] 



0.1435 
-0 . 1491 
-0.1543 



rrBq[i(U4.2I 



■0 .9972 
■0 . 9880 



fnq(kk4.3] 



0 . 9845 



Thus, from the results of the above it may be seen that a 
Multiple Correlation Coefficient better than 0.7 may be 

10 obtained using only either the features extracted from the 
power spectrum or the features extracted from the PTP; only 
part of the features offered by the method of the invention. 
Thus, it is contemplated that even higher Multiple Correla- 
tion Coefficient may be expected when these groups of fea- 

15 tures are combined in the estimation procedure. In this 

connection it should be mentioned that Multiple Correlation 
Coefficients of this size are quite satisfactory. 

EXAMPLE 3: The MORSE Esttxnator for Mixed Non- causal SARMA 
Models 

20 The concept of texture is widely recognized as being of great 
importance in many fields of application, such as, e.g., 
medical imaging, industrial quality inspection and remote 
sensing. Visual properties of a given texture are, however, 
not very tangible and many attempts has been made to extract 

25 good and meaningful features in order to describe such prop- 
erties. 

One such approach is to model the spatial correlation struc- 
ture using texture models. The two mainstream types of models 
are the Gaussian Markov Random Field (GMRF) models and the 
30 Simultaneous Auto- Regressive Moving Average (SARMA) models. 
The latter type of models is an attempt to extend the very 
powerful ARMA models (known from time- series analysis and 1-D 
signal processing) to two dimensions. 
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Whereas l-D signals are inherently causal, this is not a 
natural constraint in two dimensions. Abandoning causality, 
however, introduces severe problems in important modelling 
steps as identification and estimation. Therefore, it is very 
popular to impose an artificial directionality on data and 
thereby obtaining consistent estimates by minimizing the sum 
of squared estimation errors - i.e. Least Squares estimates. 

The analysis of non-causal SARMA models has been severely 
constrained, mainly due to the fact that the non- consistency 
of the LS estimator forces one to consider a very complicated 
Likelihood function. As a consequence, the analysis has been 
concentrated solely on pure SAR models. The toroidal ML 
estimator and an approximative ML estimator has been formu- 
lated for pure SAR models. It can be shown that the perform- 
15 ance of the approximative estimator may not be adequate when 
the considered textures is close to non-stationarity . 

Considering the exact ML estimates for toroidal SAR models 
will, in general, give very good results. It is shown, how- 
ever, that the lilcelihood function is extremely rippled in 
the vicinity of non-stationarity. This will, in many cases, 
cause the optimization routine to get 'caught' in a sub- 
optimal solution. This problem may indeed be serious as 
visual textures are often characterized by a strong positive 
correlation structure yielding estimates close to non-statio- 
25 narity. 



20 



Another problem, which is related to the Mcucimum Likelihood 
estimator and the traditional Maximum Entropy estimators, is 
that these tend to overemphasize the high-frequency content 
of the considered. 



30 Thus, a new, entropy based, estimator is introduced. Instead 
of maximizing the entropy of the associated model spectrum 
(as in MEM spectral estimation introduced in "A new algorithm 
for two-dimensional mcuciroum entropy power spectrum estima- 
tion" Lim, J. S. & Malik, N. A. (1981) IEEE Transactions on 
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Acoustics, Speech and Signal Processing, 29, 401-413), the 
entropy of the residual spectrum is mcucimized. This estimator 
is named the MORSE estimator (Mcucimum Of Residual Spectral 
Entropy) . This estimator is well behaved, even close to 
5 non-stationarity , which makes it much simpler to perform the 
non- linear optimization of the involved objective functional 
(w.r.t. the parameters) . It is possible to write the objec- 
tive functional for mixed SARMA models (i.e. including Moving 
Average parts) and obtain very good estimates. To the best of 
10 our knowledge there are no results in the literature regard- 
ing inference in non- causal SARMA models. 

The maximization of the residual spectral entropy seems, 
intuitively, to be a reasonable criterion as the underlying 
assumption of the simultaneous models is that the innovation 

15 process is white. It is obvious that the optimality of a 
given model is closely related to the chosen criterion of 
optimality, and the choice of entropy is by no means the only 
one to be imagined. Note that no assumptions are made 
concerning the driving noise (innovation process ????) 

20 sequence. We only require it to have a symmetrical 
distribution. 

The non-causal SARMA model 

Consider a wide sense stationary, zero mean bilateral SARMA 
process sampled on a square MxM grid of pixels s=(r,c)en 



ys^J2^'y»-' = ^5 * EMs-r (40) 

25 where x={r-^.r^) e N. and N. defines the order of the model. 

The only assumption that is made sibout the innovation process 
{y.} is that it is independent identically distributed random 
variables with a symmetrical density. 
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Contrary to the GMRF models, the bilateral SARMA model does 
not have to be symmetrical, but opposite neighbours must be 
associated with identical parameters (i.e. <t>^ - i>.j.. This is 
necessary to ensure the identif iability of the parameters. 
5 For the sake of simplicity, the following is restricted to 
consider only symmetrical models. 

In order to facilitate the notation when it is desired to be 
specific about the order of the model, the notation outlined 
in the following definitions is introduced: 

10 Definition 1 A rectangle with side lengths 2a+l and 2b+l, 
rotated by an angle a around the origin is denoted R(a,b,a) . 
Now the following is defined 

• R(a,b) = R{a,b, 0) 

• i?(a) = J?(a, a, 0) 

Definition 2 An ellipse with semi -major axis a and semi- 
minor axis b, rotated by an angle a around the origin is 
15 denoted E(a,b,cy). The following is defined 

• E{a,b) = E{a,b,0) 

• E(a) = E(a, a, 0) 

Thus, the order of a given model may be defined in terms of 
these geometrical shapes as, e.g., SARMA(E{1) , E (1) ) which 
means that the auto- regressive part, as well as the moving 
average part, is the four nearest neighbours. 



20 If the bilateral model is considered, the stationarity condi- 
tion is greatly weakened; all that is required is that no 
poles of the transfer function falls on the unit circle: 
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{l^ EKn^r^c" " 0} (43) 



for all J Zj- 1 =1 and |z^| = 1. 

The transfer function is given as the 2-D z-transforro of the 
impulse response function. 



and the frequency response function is obtained as the 2-D 
Fourier transform of the impulse response function 

^(/(s))=^*gf^^^-^H (45) 



where 

(46) 

^ M 2' M 2^ 



is the spatial frequency and 



cszll^r-reNj^] (47) 
) = £:Qi.[cos (27cf (s) '"r) , reW.] 



In time-series analysis and signal processing, the successful 
use of T^RMA models for spectral estimation has promoted an 
10 interest in the extension of these methods to two dimensions. 
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As for ARMA models, the 2-D parametric spectrum is propor- 
tional to the squared magnitude of the frequency response: 



S(f(s))=al\'K(f{s))\^ (48) 

where cr^^ is the variance of the innovation sequence. 

Using this approach, the problem of spectirum estimation 
5 becomes one of fitting an appropriate model to the observed 
data. Mostly, parametric spectral estimation in 2-D has been 
attempted using unilateral models. This typically leads to 
directional bias. By using a bilateral model, this particular 
problem should be solved, provided a practically applicable 
10 estimation principle is provided. Our preliminary results 
show that the MORSE estimation principle is a very useful 
tool for obtaining parametric spectral estimates without 
directional bias . 

Different approaches have been made in order to obtain esti- 
15 mates for the non-causal SAR model. There have, to the best 
of our knowledge, been no results regarding estimation of 
mixed non-causal SARMA models. Even for the pure SAR model 
the estimation of the parameters is not a simple problem. 
Below, different estimators for the non-causal SAR model will 
20 be illustrated and their shortcomings discussed. 

The Non-causal (or bilateral) Simultaneous Auto-Regressions 
(SAR) models, sometimes referred to as NCSAR, were originally 
introduced by P. Whittle in "On stationary processes in the 
plane". Biometrika, 41, 343-449. The problem of non-consist- 

25 ency of the LS estimator is addressed in this paper, and a 
large sample approach based on spectral methods was pres- 
ented. First, this method will be of use only for very simple 
models and secondly, many practical situations involve small 
san^les for which the discounting of edge effects may not be 

30 negligible. In order to illustrate how the non-consistency of 
the LS estimator affects the estimates, a simple experiment 
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is carried out. For an isotropic nearest neighbour model, 100 
realizations are generated for different parameters. In each 
case, the parameters are estimated (as if the model was 
anisotropic) . In Fig. 15, the average of the estimated para- 
5 meter is plotted versus the 'true' parameter (i.e. the para- 
- - meter with -which the texture was simulated) . It is seen that 
the LS estimator does not provide satisfactory estimates. In 
fact, the order of magnitude is off by a factor two in many 
cases . 



0 Kashyap, R. & Chellappa, R. proposed an approximation to the 
exact (toroidal) ML-estimator in "Estimation and choice of 
neighbours in spatial interaction models of image", IEEE 
Transactions on Information Theory, 29(1), 60-72. The approx- 
imation is based on a Taylor-series expansion of the deter- 

15 minant of the Jacobian (which becomes block circulant with 

circulant blocks under the torus assumption) . Only first- and 
second-order parts (in the parameters) are included. The 
improvement of the results is definitely vast, as shown in 
Fig. 16 where the above experiment is repeated. Nevertheless, 

20 as the bias is very pronounced close to non-stationarity , it 
is seen that this estimator may not be adequate in many 
practical situations where the textures tend to lie close to 
non-stationarity, and it is highly preferred to obtain a more 
suitable estimator not having these disadvantages . 



25 Grunkin investigated the exact (toroidal) ML-estimator in "On 
The Analysis of Image Data Using Simultaneous Interaction 
Models" in IMSOR, ph.d thesis # 67, Institute of Mathematical 
Statistics and Operations Research, Technical University of 
Denmark, Lyngby, 223pp. The likelihood function involved is 

30 complicated and coit^jutationally very expensive. Also, the 
likelihood function is extremely rippled in the vicinity of 
non-stationarity as shown in Fig. 17. Whether the non-linear 
optimization routine is 'caught' in a sub-optimal solution or 
not is very much dependent upon the initial guess and the 

35 parameters of the optimization algorithm (such as, e.g., step 
length) . Thus, it may take several re-estimations to be 
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convinced that the correct solution is found. This, of 
course, is merely a technical problem, and in the simple 
experiment described above, these problems can be overcome 
and the exact (toroidal) ML-estimates yield very good results 
5 as shown in Fig. 18. 

The MORSE estimator 

As mentioned above, the non- consistency of the LS estimator 
is leads to non-sensical results. The practical application 
of the exact toroidal ML estimator will, in many cases, be 
limited, as the likelihood function tends to be extremely 
rippled in the vicinity of non-stationarity . Thus, it may 
require several re- estimations to ensure that the 'true' 
optimum is found. Also, as pointed out earlier, the ML esti- 
mator may overemphasize the high- frequent content of the 
texture. In addition, the computational load involved in 
obtaining estimates for mixed SARMA models using the ML 
estimation principle is prohibitive for all practical 
purposes . 

The approach taken here is to consider the spectral content 
20 of the residuals. Finally, it will probably be extremely 
difficult to obtain ML estimates of mixed non- causal SARMA 
models . 

In this Example, a novel estimator is described that gives 
very good estimation results, has a response surface that (in 
25 our experience) is very well behaved and allows for inference 
in mixed SARMA models. 

As the underlying assumption is that the residuals should be 
white, i.e. uncorrelated, all frequencies should be represen- 
ted, with the same weight, in the spectrum of the residual 
30 process . If the residual spectrum is normalized, the entropy 
of the residual spectrum obtains its maximum when the spec- 
trum is as flat as possible. This observation is the basis of 
the Mcocimum Of Residual Spectral Entropy (MORSE) estimator. 



10 
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Consider again the difference equation defining the system: 



ys^E*ry«-r = e^^52Ms-x (49) 



The output can be viewed as the convolution of the innovation 
process with the impulse response of the system. Thus, in the 
frequency domain, the following representation is found: 

Y(.f(s)) = 'Klfis))-Eif{s)) (50) 



5 where Y(f(s)) is the Fourier transform of the output, E(f (s)) 
is the Fourier transform of the (non-observcible) innovation 
process and !K(f(s)) is the Fourier transform of the impulse 
response function, which is also known of the frequency 
response function. 

10 Using equation (50) , the power spectrum of the innovation 
process may be expressed in terms of the power spectrum of 
the observations and the parametric spectrum (assuming that 
the parameters are known). Thus it is found that: 



,^(^(.) ) JX^^L£lIliHy(f (s) ) p f^^^^^-^-^-f 
|H(£(s))|^ {i+ei^(jy^)„,,) 



(51) 



A normalization function depending upon the observations and 
15 the parameters may be defined: 

N(y,e„.e„).i;iE(/(s))p=£Jitii^.K^a| 



The parameters are subsec3[uently found by maucimizing the 
entropy of the residual spectrum. 
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Definition 3 The Maximum Of Residual Spectral Entropy 
(MORSE) estimates are found by maximizing the entropy func- 
tion lj{y.,Qpji,Qt^i^) with respect to the parameters 3.f^ and Qf^-. 



where 

(54) 



An estimate of the residual variance is obtained as 



dl=N{ir,Q.„.Q^ ) /M^ (55) 



Clearly, this function is strongly non- linear in the parame- 
ters. The optimization may be carried out using, e.g., a 
Quasi-Newton method. In this work, the IMSL implementation of 
10 the BFGS method is applied. 

Esqserimental results 

As illustrated above, the exact (torus- ) ML-estimates give 
very good results. The problem here is of a more technical 
nature. Firstly, the involved likelihood function (or objec- 
15 tive functional) is extremely complicated and computationally 
expensive to evaluate. Secondly, the likelihood function is 
very rippled in the vicinity of non-stationarity . 
Consequently any optimization-routine might easily 'get 
caught' in a sub -optimal solution. 

20 Above, importcint problems as centrality, consistency, 

efficiency of the proposed estimator have not been addressed. 
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Obviously, questions as these are difficult to answer 
analytically due to the untractable form of the involved 
expressions. Below, these cjuestions will be addressed 
empirically through simulation and estimation experiments. 

5 Estimation of pure SAR models 

In Fig. 19, the objective functional of the MORSE estimator 
is shown in the same region as the ML-estimator . It is clear 
that the objective functional of the MORSE estimator is much 
smoother - i.e. there are no ripples. It is our experience 
10 that the objective functional of the MORSE estimator never 
exhibits ripples in parameter space - not even outside the 
admissible region. This makes the non-linear optimization a 
rather straight -forward task that can be undertaken 
automatically by most standard software packages. 

15 The experiment carried out for the various estimators above 
in connection with Fig. 17 is also carried out for the MORSE 
estimator. The result is shown in Fig. 20, and it is seen 
that the results are very acceptc±>le. Note that the variance 
decreases as non-stationarity. is approached. This is due to 

20 the fact that the poles in the transfer function are 

approached. Such poles (or near-poles) create sharp peaks in 
the objective functional and thereby a much smaller variance. 

Estimation of mixed SARMA models 

A major advantage of this approach is that an estimator for 
25 mixed SARMA models is obtained 'for free' . In order to illus- 
trate that good results are obtained for the mixed SARMA 
models as well, the following experiment is conducted. A set 
of parameters are chosen for a SARMA(E{1.0) .E (l,o) ) model. 
Furthermore, these parameters are chosen not too far away 
30 from non-stationarity so that strong positive spatial auto- 
covariances are obtained. Now 100 simulations are carried out 
and the histogram of each parameter is plotted with the 
'true' parameter marked as shown in Fig. 20. 
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Another test is one where 1000 realizations of size 64x64 
pixels and 1000 realizations of size 128x128 pixels are 
generated - also not too far from non-stationarity . For each 
of these realizations, the MORSE estimator was used to obtain 
5 an estimate of the model parameters. The empirical 95%- 
confidence intervals are calculated for the mean of each 
parameter. In fig. 24, it is seen that the known values of 
the parameters fall nicely within these confidence intervals 
further indicating that the estimator is central . It should 
10 also be noted that the width of the confidence intervals 
decrease as the size of the synthesized texture samples 
increase. This indicates that the estimator is consistent. 

Thus, it is seen that very satisfactory results are obtained 
also for the mixed non-causal SARMA model. 

15 Distributional properties of the estimator 

It might intuitively be expected that the distribution of the 
estimated parameters be Gaussian if the driving noise was 
Gaussian. As this is hard to prove formally, this is made 
probable through experiments. Histograms of the parameters, 

20 estimated from the 1000 128x128 simulated textures described 
above, are shown in Fig. 25. The hypothesis that the shown 
histograms represent samples of Gaussian estimates will now 
be formally tested using the test. For each parameter, the 
estimates are divided into K classes, where K is selected as 

25 K = 1+3 . aiog^o (N) (where N is the number of observations). In 
order to make the test as strong as possible, the classes are 
selected such that the expected number of observations is 
equal in each class. This is done by generating the splits as 
the (i • 100/k) %-quantile (i=l, . . . , K-l) . The test statistic is 

30 now computed as: 




(56) 
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Where Obs^ is the number of observations in the i'th class. 
This test statistic is distributed as x^(K-3) (since both 
mean and standard deviation should be estimated) . The tests 
are summarized in the below table, where it is seen that the 
5 distributional hypothesis is easily accepted for all 
parameters . 



r. € N 


di 


mean of 0,- 


sdev of 0i 


Z 


p{z > x'm 


AR{-1,Q) 1 


-0.22314 


-0.22308 


0.00757 


5.44 


0.709 


AR{0,-1) 


-0.12118 


-0.12126 


0.00808 


6.94 


0.543 


MA{-1,0) 


0.11535 


0.11521 


0.00772 


7.512 


0.483 


MA(0,-1) 


0.24066 


0.24083 


0.00759 


6.038 


0.643 



Conclusion 

The principle on which this estimator lies is intuitively 
10 very appealing. The empirical results obtained for this 

estimator is, indeed, very promising. Using this estimator, 
also non-causal SARMA models may be considered. 

EXAMPLE 4 ; Test of the presently preferred method 

In the present example, a test was performed in order to 
15 evaluate the applicaOsility of the present method in the 

prediction of bone strength on the basis of X-ray pictures of 
the bone . 

Data set 

The data set consisted of 20 samples of the third lumbar 
20 vertebrae harvested from cadaveric specimens each of which 

were exposed to X-ray irradiation in order to obtain an. X-ray 
picture and a BHD measurement of each bone . 

Subsecj^iently , an anterior and a posterior core sample of 
trc±iecular bone was extracted from each bone, and the optimal 
25 fracture load was determined for each of these samples using 
a three-point test setup as known per s£- 
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The X-ray pictures were scanned with a spatial resolution of 
600 lines per cm. In the final resolution, an area is 
selected centrally positioned between the positions of the 
drilled-out samples. The positions of the core sanples may be 
5 partly or fully comprised in the selected area. 

It was decided to choose the mean optimal fracture load of 
the two core samples of a bone as the strength of the bone. 
This largely corresponds to the strength of the bone in the 
central part as determined by linear interpolation of the 
10 strengths of the two samples. The results of the present test 
seems to confirm that this decision is reasonable. 

-P-re-Drocessina of the imacre material 

In order to have as good images as possible in the final, 
much smaller resolution, the scanned X-ray images had a large 
15 spatial resolution (600 lines per centimetre) . 

The present images were so large that a scale reduction using 
a Burt-Adelson Gaussian pyramid was impossible due to the 
computer having too little memory. Thus,, as a preliminary 
reduction was performed by firstly smoothing the image using 
20 a 3x3 mean filter and subsequently discarding every other 
pixel in the image . 

After this preliminary reduction, the image had a size 
reducible by the Burt-Adelson Gaussian pyramid to obtain the 
final resolution (150 lines per centimetre) . 

25 Subsequent to the reduction of the image, the image was 

background corrected by subtracting a 31x31 median filtered 
version of itself . 

From the corrected image, a Region Of Interest (ROD having a 
size of 256x256 pixels was selected- PrefercQjly, this ROI was 
30 selected close to the center of the bone due to the fact that 
a large gradient of the thickness of the trabeculae was 
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present close to the anterior and posterior ends of the bone. 
This effect and the fact that the bone is rounded towards its 
ends created a number of non- linear effects in the image 
dynamics which cannot be dealt with in a trivial manner. 

5 Having extracted the Region Of Interest, the image "was 
filtered in order to remove high frequency noise. As the 
signal-to-noise relation of the present image was very high, 
the high frequency noise may be removed in a "gentle" manner: 

fourtr > fft.$$ 
10 htof -CR M < fft.$$ I powerpix -p 1.8 > mag.-$$ 
htof -CR P < fft.$$ > pha.$$ 

combine -pm pha.$$ < mag.$$ j inv. fourtr | htof -CR R j scale 
rm * .$$ 

This manner is an alternative to the method of algorithm l on 
15 pages 27 and 28, as, in this method, T=2 and p=1.8. 

From this pre-processed image, the Projected Tredaecular 
Pattern (PTP) is extracted in the manner described above in 
connection with Algorithm 2. Substantially all features 
described in Example 1 were calculated. 

20 Statistical methods 

It was attempted to construct a General Linear Model for 
predicting the ultimate strength. In this connection, it is 
imperative to identify the inportant parameters and eliminate 
the unimportant ones . As there are a large number of 

25 variables and as only a small number of observations were 
evaluated, it was not possible to directly use a Backwards 
Elimination method as (X^X) typically will be singular in 
this situation. In addition, there is also the risk of 
eliminating the wrong variables. Therefore, in the following 

3 0 a method is used which resembles a Forward Elimination method 
performed as follows: 
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1. Search through the list of possible features until the 
first significant feature is identified. 

2. Continue the search until another significant feature is 
identified and included in the model. 

5 3. Test whether previous features should be eliminated. 

Steps 2 . cuid 3 . are repeated until all features have been 
investigated. Naturally, a level of significance should be 
decided for the described tests . In the present 
investigation, this level is chosen at p=0.10, as the amount 
10 of data does not make a very conservative test possible. 

It became clear that the optimal fracture load related 
exponentially to the parameters of the model, whereby the 
model was performed as a regression with the logarithm to the 
optimal fracture load as a dependent variable. 

15 Choosing the parameters according to the above principle, 
only three significant parameters were found: § (See 
Algorithm 4) , Formula (23) squared, and fj.^^ (See Formula 
(3)). 

In the following, these parameters are evaluated in two 
20 models, one using age as a variable and one not using age. 

Afodellinq without acre 

A General Linear Model not using age as a descriptive 
variable was made. This model will, thus, seek to predict the 
optimal fracture load on the basis of the extracted 
25 structural indices . The result is summarized in the below 
program printout : 

Residuals : 

Min IQ Medicui 3Q Meoc 

-0.4327 -0.2409 0.03446 0.1625 0.5578 



) 
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(Formula 23 )_ 



Value 

11.2930 
-70 . 8473 
-0 . 0287 
143625 .7515 



Std. Error t value Pr(>|t|) 
1 . 8815 
21.1912 
0.0037 
-53274.0853 



6.0020 0.0000 

3.3432 0.0044 

■7.8621 0.0000 

2.6960 0.0166 



10 



Residual standard error: 0,2901 on 15 degrees of freedom 
Multiple R-Squared: 0.8158 

F- statistic: 22.14 on 3 and 15 degrees of freedom, the p- 
value is 9.17e-06 



Correlation of Coefficients 

(Intercept) 
a -0.9873 



Ms 



15 (Formula 23) ^ 

Explaining model : 



-0.5412 
0 .4347 



0 .4485 
-0 .5455 



Ms' 



-0 . 1318 



F^=exp { -7 0 . 847 3 • (Formuia23 ) 2-0 . 0287 -5 
♦143625 .7 515 -nl+ll. 2930) 



(57) 



|No. 


log (Msrd) 


log (Predicted) 


Residual | 


238 


4.574778 


4.302492 


0.272286112 


240 


4.259244 


4.168912 


0.090332350 


242 


4.516842 


4.742086 


-0.225244419 


247 


4.677491 


4.800312 


.0.122820941 


250 


3.508197 


3.473737 


0.034458077 


254 


4.070737 


4.402031 


.0.331293705 


258 


2.601764 


2.864725 


-0.362961146 


259 


3.724428 


3.682320 


0.042107301 


265 


4.443404 


3.885623 


0.557781393 


267 


5.382806 


5.098572 


0.284233332 


269 


4.099266 


3.949321 


0.149945096 


279 


4.434097 


4.125558 


0.308539021 


280 


3.844675 


4.277384 


.0.432709063 


281 


3.372283 


3.310585 


0.061698025 


289 


3.823039 


3.808975 


0.014063294 


293 


3.779999 


3.604959 


0.175040294 


296 


3.527713 


3.784313 


-0.256600239 


298 


3.896777 


3.886981 


0.009796051 


313 


4.397531 


4.666183 


.0.268651834 
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In Fig. 26 a plot of the logarithmic optimal fracture load 
against the predicted strength is illustrated. In this 
figure, the measurement points are represented by their 
number in the datod3ase . 

5 At this point it is, naturally, interesting to know whether 
the BMD measurements performed on all samples may take part 
as a statistically significant feature contributing to 
increasing the multiple R- squared. The result may be seen 
from the below printout 



10 Residuals 
Min 

-0 .4342 



IQ 

-0.2389 



Median 
0 . 0385 



3Q 

0 ,1615 



Max 

0 .5598 



Coefficients 
15 (Intercept) 

a 

Formula 23 
squared 
20 BMD 



Value 

11.3316 
-71.1532 

-0.0288 
144161 . 1908 

-0 .0190 



Std. Error t value Pr(>|t|) 

2.2940 4.9397 0.0002 

23.9423 -2.9719 O.OlOl 

0.0044 -6.5255 0.0000 

57644.4790 2.5009 0.0254 

0.5953 -0.0319 0.9750 



Residual standard error: 0.3002 on 14 degrees of freedom 
Multiple R-Squared: 0.8158 

F-statistic: 15.5 on 4 and 14 degrees of freedom, the p-value 
is 8.43e-05 



25 Correlation of Coefficients: 

(Intercept) S 



Ms 

30 Formula 23 
squared 
BMD 



-0.9798 
-0.6661 
0 .5070 

-0.5285 



0 .5587 
-0 .5949 

0.4009 



Ms 



Formula 23 
squared 



-0.2583 



0.5154 -0.2914 
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It is obvious that BMD is not a significant feature and that 
it can be eliminated from the model. Using BMD as the only 
descriptive parameter, a multiple R-squared of R^=0.2142 is 
obtained indicating that BMD is a poor indicator of the 
5 quality of trabecular bone. Thus, the aspects of fracture 

risk, relating to trabecular bone are poorly described by BMD. 



Modelling with aae 



Using age as a describing parameter, the following model is 
obtained : 



10 Residuals: 

Min IQ Median 3Q Max 

-0.3777 -0.1385 0.02544 0.1517 0.328 



Coefficients: 

Value Std. Error t value Pr(>|t|) 

15 (Intercept) 11.8882 1.6616 7.1545 0.0000 

a -68.3085 18.5326 -3.6858 0.0024 

Ms^ -0.0236 0.0039 -6.1152 0.0000 

Formula 23 152951.1278 46677.8586 3.2767 0.0055 
squared 

20 AGE -0.0160 0.0067 -2.3827 0.0319 



Residual standard error: 0.2533 on 14 degrees of freedom 
Multiple R-Squared: 0.8689 

F-statistic: 23.2 on 4 and 14 degrees of freedom, the p-value 
is 4.707e-06 

Formula 23 
squared 
f 

-0.0617 

-0.5613 -0.0838 



25 Correlation of Coefficients: 

(Intercept) oL 

§ -0.9658 

-0.3585 0.4028 
30 Formula 23 0.4408 -0.5378 
squared 

AGE -0.1503 -0.0575 
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Explaining model : 



f^^=exp (-68 . 3 085- (ForintjIa23) 2-0 . 0236-5 

+152951. 127 8-ji|-0 .016 rAGSf 11. 8882) 



No. 


log(Msrd) 


log (Predicted) 


Residual 


238 


4.574778 


4.293879 


0.280898703 


240 


4.259244 


4.262127 


-0.002882869 


242 


4.516842 


4.505306 


0.011535494 


247 


4.677491 


4.766867 


•0.089375665 


250 


3.508197 


3.413527 


0.094669315 


254 


4.070737 


4.435918 


•0.365180333 


258 


2.501764 


2.852028 


•0.350263751 


259 


3.724428 


3.625144 


0.099283320 


265 


4.443404 


4.238432 


0.204971765 


267 


5.382806 


5.273592 


0.109213882 


269 


4.099266 


3.964857 


0.134408820 


279 


4.434097 


4.106100 


0.327997388 


280 


3.844675 


4.032239 


•0.187563938 


281 


3.372283 


3.394572 


•0.022288673 


289 


3.823039 


3-797599 


0.025439956 


293 


3.779999 


3.610914 


0.169084738 


296 


3.527713 


3.829923 


-0.302210802 


298 


3.896777 


3.656788 


0.239989248 


313 


4.397531 


4.775257 


-0.377726598 



In Fig. 27 a plot of the logarithmic optimal fracture load 
against the predicted strength is illustrated. In this 
5 figure, the measurement points are represented by their 
number in the database. 

Discussion 

From the result of the present example, it is clear that the 
optimal fracture load of trabecular bone may be predicted 

10 using only textural parameters - and even with a high 

precision. Three image features are used: Formula 23 squared, 
/ig^/ aj^d, a describing the density of the horizontal 
trsdseculae (the thickest) , the break down of the trabecular 
structure, and the global treOaecular density, respectively. 

15 In addition, it is seen that neither sex nor BMD can add 
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information to the description obtained using the cibove three 
image features . 

However, it is seen that age is a significant parameter. 
Including age in the model increases the multiple R- squared 
5 to R^=0.8689.- It should be noted that the two extreme 

observations (258 and 267) contribute quite a lot to the 
value of R^. Omitting these observations gives R^=0.7 53B, 
including age, which should still be regarded as a high 
multiple R- squared. 

10 However, as is described in connection with Algorithm 1, the 
parameters used in restoring the image will a large effect on 
the applicability of the features extracted from the 
resulting image. In fact, using p=1.4 and T=4 instead of 
p=T=2 in algorithm l, the feature a^^ alone gives a R^ of 

15 0.92. 

In Fig. 28, the correlation between and the logarithm to 
the optimal fracture load is illustrated for a test based on 
the picture and optimal fracture load - material described 
above. Even if the extreme point in the lower right corner of 
20 the figure is removed, a correlation of 0.88 is obtained. 

Thus, a very large difference may be seen depending on the 
actual manner of obtaining the PTP of the image. 
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1. A method for estimating the bone quality of a vertebrate, 
on the basis of two-dimensional image data comprising infor- 
mation relating to the treibecular structure of at least a 
5 part of a bone of the vertebrate, the image data being data 
obtained by exposing at least the part of the bone to elec- 
tromagnetic radiation, the method comprising sxibjecting the 
image data to a statistical analysis comprising 

a background correction procedure in which low frequent 
0 intensity variations not related to the trabecular 

structure of the bone is reduced relative to image data 
related to the trabecular structure of the part of the 
bone, 

an image manipulation and feature extraction procedure 
5 wherein at least the local image intensity information 

as well as variation in the local intensity are utilized 
to extract information related to the treibecular 
structure of the part of the bone, the image 
manipulation and feature extraction procedure comprising 
0 subjecting the image data to at least one of the 

following procedures: 

(a) obtaining an estimate of the parametric estimate 
of the power spectrum of the image data and 
extracting features relating to the energy 

5 distribution of the parametric estimate, 

(b) obtaining, on the basis of image data on which a 
Fourier method has been used to emphasize the 
information in the image data relating to the 
trabecular structure, an estimate of a grey- level 

0 co-occurrence matrix and extracting at least one 

feature on the basis of the estimated co- 
occurrence matrix, 

(c) obtaining an estimate of the projected trabecular 
pattern of the image data by using a Fourier 

5 method to eir^hasize the information in the image 

data relating to the traibecular structure and 
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slibjecting the manipulated image data to a 
morphological operation, and extracting features 
relating to the trabecular structure from the 
estimated projected trabecular pattern, 



5 



(d) obtaining, on the basis of a frequency analysis 



of the image data, features relating to the 
periodicity of the trabecular strxicture of the 
part of the bone. 



10 



and an estimation procedure in which the bone quality of 
the vertebrate is estimated on the basis of the derived 
features and optionally other features related to the 
bone or the vertebrate and a predetermined relationship 
between the features and reference bone quality parame- 



ters . 



15 2. A method according to claim 1, wherein the electromagnetic 
radiation is X-ray radiation. 

3 . A method according to any of the preceding claims , wherein 
the image data is scanned from an X-ray film. 

4. - A method according to claim 3, wherein the resolution of 
20 the X-ray film is at least 4 pairs of lines per centimetre, 

such as at least 5 pairs of lines per centimetre. 

5. A method according to claim 4, wherein the resolution of 
the X-ray film is at least 10 pairs of lines per centimetre, 
such as at least 25 pairs of lines per centimetre, preferably 

25 at least 50 pairs- of lines per centimetre. 

6 . A method according to claim 5 , wherein the resolution of 
the X-ray film is at least 100 pairs of lines per centimetre, 
such as at least 250 pairs of lines per centimetre, preferoib- 
ly at least 500 pairs of lines per centimetre, such as at 

30 least 600 pairs of lines per centimetre. 

7 . A method according to claim 3 , wherein t.he scanning has 
been performed at a resolution of at least 10 lines per cm. 
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8. A method according to claim 7, wherein the resolution is 
at least 25 lines per cm. 

9. A method according to claim 8, wherein the resolution is 
at least 100 lines per cm, such as at least 200 lines per cm, 

5 such as at least 250 lines per cm. 

10. A method according to any of the preceding claims, where- 
in the background correction procedure comprises at least 
reducing or optionally removing low frequency information 
having a frequency significantly lower than the spacing of 

10 the projected trabeculae. 

11. A method according to claim 10, wherein information 
having frequencies half or less than the spacing of the 
projected trabeculae is at least reduced or optionally 
removed . 

15 12. A method according to claim 11, wherein information 

having frequencies being a quarter or less than the spacing 
of the projected trabeculae is at least reduced or optionally 
removed . 

13. A method according to claim 12, wherein information 

2 0 having frequencies being a tenth or less the spacing of the 
projected trabeculae is at least reduced or optionally 
removed . 

14. A method according to any of claims 10-13, wherein the 
background correction procedure comprises generating second- 

25 ary image data as a result of performing a medism filtering 
with a predetermined kernel size and subtracting this result 
from the original image data. 

15. A method according to any of claims 10-13, wherein the 
background correction procedure comprises generating second- 

30 ary image data as a result of performing a mean filtering 
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with a predetermined kernel size and subtracting this result 
from the original image data. 

16. A method according to any of claims 10-13, wherein the 
background correction procedure comprises globally fitting a 

5 two-dimensional polynomial to the image data and generating 
background corrected image data on the basis of the residuals 
of the fitting procedure. 

17. A method according to claim 16, wherein the order of the 
polynomial is at the most 15, such as at the most 10, prefer- ' 

10 ably at the most 5. 

18. A method according to claim 14 or 15, wherein the kernel 
size is at the most 1/2 of the image data, such as at the 
most 1/4 of the image data, preferably at the most l/lO of 
the image data, more preferably at the most 1/20 of the image 

15 data. 

19. A method according to any of the preceding claims, where- 
in the parametric estimate of the power spectrum is estimated 
by using a method selected from the group consisting of 
direct methods, auto-covariance methods and parametric 

20 methods. 

20. A method according to claim 19, wherein the image data, 
optionally weighted with a window, is subjected to a Fast 
Fourier Transformation in order to generate the parametric 
estimate . 

25 21. A method according to claim 19, wherein the auto- 
covariance function estimated on the basis of the image data, 
optionally weighted with a window, is subjected to a Fast 
Fourier Transformation in order to generate the parametric 
estimate . 

30 22 . A method according to claim 19, wherein one of the 
methods on the basis of which the parametric estimate is 
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Obtained is chosen from the group consisting of: causal 
Simultaneous Auto-Regressive Moving Average (SARMA) models, 
non-causal SARMA models, Gaussian Markov Random Field models 
and Maximum Entropy Spectral Estimates and wherein the order 
5 of the model is identified, the parameters of the model are 
estimated and the spectrum of the estimated model is gener- 
ated. 

23. A method according to claim 22, wherein the method on the 
basis of which the parametric estimate is obtained, is a non- 
10 causal SARMA model using an estimator according to claim 50. 

24. A method according to any of the preceding claims, where- 
in at least one feature is related to parameters of a contour 
encompassing at least a predetermined percentage of the 
energy of the power spectrum. 

25. A method according to claim 24, wherein the contour is 
determined by 

(a) defining a contour, in terms of one axes in each 
dimension, around the centre of the power spectrum, all 
points on the contour having substantially the same 
distance to the centre of the power spectrum, and the 
contour encompassing less than the predetermined 
percentage of the energy of the power spectrum, 

(b) for each dimension of the data calculating the per- 
centage of the energy of the power spectrum encompassed 
by a dilated contour in which the axis in the dimension 
in question is increased by a predetermined distance, 

(c) increasing the axis of the contour, in the dimension in 
which the increase in energy encompassed by the dilated 
contour is the largest, by the predetermined distance 
and 

(d) repeating steps (b) euid (c) until the percentage 
encompassed by the contour exceeds or equals the 
predetermined percentage . 



15 



20 



25 
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26. A method according to claim 24 or 25, wherein the contour 
is ellipsoidal, optionally spherical. 

27. A method according to claim 25 or 26, wherein the axis of 
the contour are orthogonal and are defined along principal 

5 directions of the power spectrum. 

28. A method according to claim 1, wherein procedure (c) is 
used and wherein the Fourier method comprises 

subjecting the image data to Fourier transformation, 
performing a subsequent mathematical transformation of 
10 the Fourier transformed data, 

converting the transformed data back into the spatial 
domain . 

29. A method according to claim 28, wherein the subsequent 
mathematical transformation is linear or non- linear. 

15 30. A method according to claim 29, wherein the subsequent 
mathematical transformation is performed by raising the data 
to a power larger than 1, such as a power in the range of 
1.1-10, preferably 1.5-4, such as aibout 2.0. 

31. A method according to any of claims 28-30, wherein a 
20 third mathematical transformation is performed in order to 

preserve only part of the magnitude information of the trans- 
formed data. 

32. A method according to any of claims 28-31, wherein the 
emphasized information in the image data relating to the 

25 trabecular structure is subjected to a grey- scale morphologi- 
cal operation. 

33. A method according to claim 32, wherein the grey-scale 
morphological operation is a top-hat operation wherein all 
pixels fitting in the top-hat are assigned a value different 

30 from a specified background pixel value. 
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34. A method according to claim 33, wherein the image data 
which has been subjected to the top-hat operation is thres- 
holded so as to assign a foreground pixel value to the pixels 
having a pixel value different from the background pixel 

5 value . 

35. A method according to any of claims 32-34, wherein the 
image data which has been subjected to the grey- scale mor- 
phological operation is subjected to further noise reduction 
to substantially remove isolated single pixels and pixel 

10 groups smaller than a predetermined threshold number of 
pixels having a pixel value different from the background 
pixel value - 

36. A method according to any of claims 32-35, wherein an 
operation is performed comprising determining, for each 

15 background pixel, the distance according to a given metric 
from the background pixel to the nearest foreground pixel, 
and wherein at least one feature is extracted from the result 
of the operation. 

37. A method according to claim 36, wherein at least one 
20 feature relates to a mean value and/or standard deviation 

and/or the coefficient of variation and/or skewness and/or 
kurtosis of the determined distances from all background 
pixels to the nearest foreground pixel. 

38. A method according to any of claims 32-35, wherein an 
25 operation is performed cotrprising determining a measure for 

each background pixel based upon determining the distance 
from the background pixel in a number of given directions in 
the image data to the nearest foreground pixel and at least 
one feature is extracted from the result of the operation. 

30 39. A method according to claim 38, wherein at least one 
feature relates to a mecui value and/or standard deviation 
and/or the coefficient of variation and/or skewness and/or 
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kurtosis of the determined measures for all background 
pixels . 

40. A method according to claim 38 or 39, wherein the measure 
for each background pixel is a mean value and/or standard 
deviation and/or the coefficient of variation and/or skewness 
and/or kurtosis and/or maximum distance of the determined 
distances in the given directions. 

41. A method according to any of the preceding claims, 
wherein at least one feature related to the bone or to the 

) vertebrate, such as age and/or sex, and/or species, and/or 
race and/or the specific bone considered in the vertebrate, 
and/or an estimated Bone Mineral Density of the bone, and/or 
an estimated Bone Mineral Content of the bone, is included in 
the estimation procedure. 

15 42. A method according to claim 41. wherein the Bone Mineral 
Density is estimated by including data from a reference 
object in the exposure of the bone to the electromagnetic 
radiation and where the Bone Mineral Density is estimated by 
correlating the absorption of the electromagnetic radiation 

20 of the bone and of the reference object. 

43. A method according to any of the preceding claims, 
wherein the estimation procedure is based on a statistical 
model, taking into account the correlation structure in the 
data set, so as to assign appropriate weights to the 

25 significant features in accordance with the predetermined 
relationship . 

44. A method according to any of the preceding claims, 
wherein the reference bone quality is an absolute or relative 
bone quality of bones for which image data have been 

30 processed in accordance with any of the preceding claims. 

45. A method according to claim 44, wherein the reference 
bone quality is determined on the basis of measurement of the 
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mechanical bone strength, and/or Bone Mineral Density 
measurement, and/or on the basis of a Bone Mineral Content 
measurement, and/or on the basis of a score of the bone 
strength by a skilled radiologist . 

46. A method according to any of the preceding claims, 
wherein the predetermined relationship between the features 
and reference bone quality parameters is established on the 
basis of reference bone quality parameters and features 
extracted according to any of the preceding claims. 

47. A method according to claim 46, wherein the predetermined 
relationship is defined in terms of a model taken from the 
group consisting of: a General Linear Model, a Generalized 
Linear Model, an Artificial Neural Network, a Causal 
Probabilistic Net and Classification And Regression Trees. 

15 48. A method according to any of the preceding claims, 
wherein the vertebrate is a human. 

49. A method according to any of the preceding claims, 
wherein the vertebrate is a horse. 



10 



20 



50. A method according to any of the preceding claims, 
wherein the vertebrate is a large ape, a great ape or an 
anthropoid ape . 

51. A method according to any of the preceding claims, 
wherein the vertebrate is a pig. 

52. A method according to any of the preceding claims, 
25 wherein the vertebrate is a cow. 



53. A method according to any of the preceding claims, 
wherein the bone is taken from the group consisting of 
radius, femur, corpus vertebrae (LI, L2 , L3 , L4 , L5, Tl , T2, 
T3, T4. T5, T6. T7 , T8 , T9, TIO, Til, T12. CI, C2 , C3 , C4 , 
30 C5, C6, C7) , calcaneus, talus, os carpi, metatars, metacarpi. 
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falanges, tibia, fibula, patella, ulna, humerus, mandible, 
clavicula, scapula, os coxae, os naviculare, os cuboideum, os 
cuneiform I, os cuneiform II, os cuneiform III, os sacrum, os 
coccygis . 

-5 54 . A method for estimating the bone quality of a vertebrate, 
on the basis of two-dimensional image data comprising 
information relating to the trabecular structure of at least 
a part of a bone of the vertebrate, the image data being data 
obtained by exposing at least the part of the bone to 

10 electromagnetic radiation, the method comprising subjecting 
the image data to a statistical analysis comprising 



15 



20 



a background correction procedure in which low frequent 
intensity variations not related to the trsQDecular 
structure of the bone is reduced relative to image data 
related to the trabecular structure of the part, of the 
bone , 

an image manipulation and feature extraction procedure 
comprising subjecting the image data to at least one of 
the following procedures : 

(a) obtaining an estimate of the parametric estimate 



of the power spectrum of the image data and 
extracting features relating to the energy 
distribution of the parametric estimate , 



(b) obtaining, on the basis of image data on which a 



25 



30 



Fourier method has been used to emphasize the 
information in the image data relating to the 
trcibecular structure, an estimate of a grey-level 
co-occurrence matrix and extracting at least one 
feature on the basis of the estimated co- 
occurrence matrix. 



(c) obtaining an estimate of the projected trabecular 



35 



pattern of the image data by using a Fourier 
method to emphasize the information in the image 
data relating to the trabecular structure and 
subjecting the manipulated image data to a 
morphological operation, and extracting features 
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relating to the trabecular structure from the 

estimated projected trabecular pattern, 
(d) obtaining, on the basis of a frequency analysis 

of the image data, features relating to the 
5 periodicity of the trabecular strxicture of the 

part of the bone, 
and an estimation procedure in which the bone quality of 
the vertebrate is estimated with a Multiple Correlation 
Coefficient better than 0.5 on the basis of the derived 
10 features and optionally other features related to the 

bone or the vertebrate and a predetermined relationship 
between the features and reference bone quality 
parameters . 



55. A method according to claim 54, wherein the bone quality 
15 of the vertebrate is estimated with a Multiple Correlation 

Coefficient better than 0.55, such as better than 0.6, 
preferedaly better than 0.65, such as better than 0.7, 
prefereUaly better than 0.8, such as better than 0.85. 

56. A method for estimating Non-Causal Simultaneous Auto- 
20 Regressive Moving Average models in two or more dimensions, 

the method comprising optimizing a given direct measure of 
the flatness of the residual spectrum of the model . 

57. A method according to claim 56, wherein the optimization 
of the given flatness measure is obtained by 

25 (a) generating a set of initial parameters for the model, 

(b) generating the residual spectrum of the model on the 
basis of the parameters, 

(c) obtaining the measure of the flatness of the residual 
spectrum, 

3 0 (d) obtaining a new iterate of the parameters on the basis 
of the flatness measure and a search direction in 
parameter space, 
(e) repeating steps (b)-{d) until given stop criterion is 
reached . 
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58. A method according to claims 56 or 57, wherein the 
measure of flatness of the residual spectrum is obtained as 
the entropy of the normalized residual power spectrum, dis- 
regarding the DC value . 
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AMENDED CLAIMS 

[received by the International Bureau on 13 February 1996 (13.02.96); 

original claims 1-58 replaced by 
amended claims 1-60 (12 pages)] 
1. A method for estimating the bone quality of a vertebrate, 
on the basis of two-dimensional image data comprising infor- 
mation relating to the trabecular structure of at least a 
part of a bone of the vertebrate, the image data being data 
obtained by exposing at least the part of the bone to elec- 
tromagnetic radiation, the method con^jrising subjecting the 
image data to a statistical analysis comprising 

a background correction procedure in which low frequent 
intensity variations not related to the trabecular struc- 
ture of the bone is reduced relative to image data 
related to the trabecular structure of the part of the 
bone, 

an image manipulation and feature extraction procedure 
wherein at least the local image intensity information as 
well as variation in the local intensity are utilized to 
extract information related to the trabecular structure 
of the part of the bone, the image manipulation and fea- 
ture extraction procedure comprising subjecting the image 
data to at least one of the following procedures: 
(a) obtaining an estimate of the parametric estimate 



5 



0 



5 



(b) 



(c) 



of the power spectrum of the image data and 
extracting features relating to the energy 
distribution of the parametric estimate, 
obtaining, on the basis of image data on which a 
Fourier method has been used to emphasize the 
information in the image data relating to the 
trabecular structure, an estimate of a grey- level 
co-occurrence matrix and extracting at least one 
feature on the basis of the estimated co- 
occurrence matrix, 

obtaining an estimate of the projected trabecular 
pattern of the image data by using a Fourier 
method to emphasize the information in the image 
data relating to the trabecular structure and 
subjecting the manipulated image data to a 
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morphological operation, and extracting features 
relating to the treUaecular structure from the 
estimated projected trabecular pattern. 



(d) obtaining, on the basis of a frequency analysis 



10 



5 



of the power spectrum of the image data, a 
smoothed periodogram of the image data, or an 
estimate of the parametric estimate of the power 
spectrum of the image data, features relating to 
the periodicity of the trabecular structure of 
the part of the bone. 



15 



and an estimation procedure in which the bone quality of 
the vertebrate is estimated on the basis of the derived 
features and optionally other features related to the 
bone or the vertebrate and a predetermined relationship 
between the features and reference bone quality parame- 



ters . 



2. A method according to claim l, wherein the electromagnetic 
radiation is X-ray radiation. 

3. A method according to any of the preceding claims, wherein 
20 the image data is scanned from an X-ray film. 

4 . A method according to claim 3 , wherein the resolution of 
the X-ray film is at least 4 pairs of lines per centimetre, 
such as at least 5 pairs of lines per centimetre, 

5. A method according to claim 4, wherein the resolution of 
25 the X-ray film is at least 10 pairs of lines per centimetre, 

such as at least 25 pairs of lines per centimetre, preferably 
at least 50 pairs of lines per centimetre. 

6. A method according to claim 5, wherein the resolution of 
the X-ray film is at least 100 pairs of lines per centimetre, 

30 such as at least 250 pairs of lines per centimetre, prefereJa- 
ly at least 500 pairs of lines per centimetre, such as at 
least 600 pairs of lines par centimetre. 
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7. A method according to claim 3, wherein the scanning has 
been performed at a resolution of at least 10 lines per cm. 

8. A method according to claim 7, wherein the resolution is 
at least 25 lines per cm. 

5 9. A method according to claim 8, wherein the resolution is 
at least 100 lines per cm, such as at least 200 lines per cm, 
such as at least 250 lines per cm. 

10. A method according to any of the preceding claims, where- 
in the background correction procedure comprises at least 

10 reducing or optionally removing low frequency information 
having a frequency significantly lower than the spacing of 
the projected trabeculae. 

11. A method according to claim lo, wherein information 
having frequencies half or less than the spacing of the 

15 projected trabeculae is at least reduced or optionally 
removed . 



12. A method according to claim 11, wherein information 
having frequencies being a quarter or less than the spacing 
of the projected trabeculae is at least reduced or optionally 

2 0 removed . 

13. A method according to claim 12, wherein information 
having frequencies being a tenth or less the spacing of the 
projected trc±>eculae is at least reduced or optionally 
removed . 



25 14. A method according to any of claims 10-13, wherein the 
background orrection procedure comprises generating second- 
ary image d .-.a as a result of performing a median filtering 
with a pr .• cermined kernel size and subtracting this result 
from the ginal image data. 
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15. A method according to any of claims 10-13, wherein the 
background correction procedure comprises generating second- 
ary image data as a result of performing a mean filtering 
with a predetermined kernel size and subtracting this result 

5 from the original image data. 

16. A method according to any of claims 10-13, wherein the 
background correction procedure con5)rises globally fitting a 
two-dimensional polynomial to the image data and generating 
background corrected image data on the basis of the residuals 

10 of the fitting procedure. 

17. A method according to claim 16, wherein the order of the 
polynomial is at the most 15, such as at the most 10, prefer- 
ably at the most 5 . 

18. A method according to claim 14 or 15, wherein the kernel 
15 size is at the most 1/2 of the image data, such as at the 

most 1/4 of the image data, preferably at the most 1/10 of 
the image data, more preferably at the most 1/20 of the image 
data. 

19. A method according to any of the preceding claims, where - 
20 in the parametric estimate of the power spectrum is estimated 

by using a method selected from the group consisting of 
direct methods, auto- covariance methods and parcunetric 
methods . 

20. A method according to claim 19, wherein the image data, 
25 optionally weighted with a window, is subjected to a Fast 

Fourier Transformation in order to generate the parametric 
estimate. 

21. A method according to claim 19, wherein the auto- 
covariance function estimated on the basis of the image data, 

30 optionally weighted with a window, is subjected to a Fast 
Fourier Transformation in order to generate the parametric 
estimate. 
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22 . A method according to claim 19 , wherein one of the 
methods on the basis of which the paraunetric estimate is 
obtained is chosen from the group consisting of: causal 
Simultaneous Auto -Regressive Moving Average (SARMA) models, 
5 non- causal SARMA models, Gaussian Markov Random Field models 
and Mcucimum Entropy Spectral Estimates and wherein the order 
of the model is identified, the parameters of the model are 
estimated and the spectrum of the estimated model is gener- 
ated. 

10 23. A method according to claim 22, wherein the method on the 
basis of which the parametric estimate is obtained, is a non- 
causal SARMA model using an estimator according to claim 58. 

24. A method according to any of the preceding claims, where- 
in at least one feature is related to parameters of a contour 

15 encompassing at least a predetermined percentage of the 
energy of the power spectrum. 

25. A method according to claim 24, wherein the contour is 
determined by 

(a) defining a contour, in terms of one axes in each 

20 dimension, around the centre of the power spectrum. 



25 



30 



(b) 



(c) 



all points on the contour having substantially the 
same distance to the centre of the power spectrum, 
and the contour encompassing less than the predeter- 
mined percentage of the energy of the power spectrum, 
for each dimension of the data calculating the per- 
centage of the energy of the power spectrum 
encompassed by a dilated contour in which the cucis in 
the dimension in question is increased by a predeter- 
mined distance, 

increasing the axis of the contour, in the dimension 
in which the increase in energy encompassed by the 
dilated contour is the largest, by the predetermined 
distance and 
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(d) repeating steps (b) and (c) until the percentage 

encompassed by the contour exceeds or equals the 
predetermined percentage . 

26. A method according to claim 24 or 25, wherein the contour 
5 is ellipsoidal, optionally spherical! 

27. A method according to claim 25 or 26, wherein the axis of 
the contour are orthogonal and are defined along principal 
directions of the power spectrum. 



28. A method according to claim l, wherein procedure (c) is 
0 used and wherein the Fourier method comprises 

subjecting the image data to Fourier transformation, 
performing a subsequent mathematical transformation of 
the Fourier transformed data, 

converting the transfoirmed data back into the spatial 
5 domain . 

29. A method according to claim 28, wherein the sxabsequent 
mathematical transformation is linear or non- linear. 

30. A method according to claim 29, wherein the sxabsequent 
mathematical transformation is performed by raising the data 

0 to a power larger than 1, such as a power in the range of 
1.1-10, preferably 1.5-4, such as about 2.0. 

31. A method according to any of claims 28-30, wherein a 
third mathematical transformation is performed in order to 
preserve only part of the magnitude information of the trans - 

5 formed data. 

32. A method according to any of claims 28-31, wherein the 
emphasized information in the image data relating to the 
trabecular stjructure is subjected to a grey- scale morphologi- 
cal operation. 
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33. A method according to claim 32, wherein the grey-scale 
morphological operation is a top-hat operation wherein all 
pixels fitting in the top-hat are assigned a value different 
from a specified background pixel value. 

5 34 . A method according to claim 33, wherein the image data 
which has been subjected to the top -hat operation is thres- 
holded so as to assign a foreground pixel value to the pixels 
having a pixel value different from the background pixel 
value . 

35. A method according to any of claims 32-34, wherein the 
image data which has been siibjected to the grey- scale mor- 
phological operation is subjected to further noise reduction 
to substantially remove isolated single pixels and pixel 
groups smaller than a predetermined threshold number of 
pixels having a pixel value different from the background 
pixel value. 

36. A method according to ajiy of claims 32-35, wherein an 
operation is performed comprising determining, for each 
background pixel, the distance according to a given metric 

20 from the background pixel to the nearest foreground pixel, 

and wherein at least one feature is extracted from the result 
of the operation. 

37. A method according to claim 36, wherein at least one 
feature relates to a mean value and/or standard deviation 

25 and/or the coefficient of variation and/or skewness and/or 
kurtosis of the determined distances from all background 
pixels to the nearest foreground pixel. 

38. A method according to any of claims 32-35, wherein an 
operation is performed comprising determining a measure for 

30 each background pixel based upon determining the distance 

from the background pixel in a nximber of given directions in 
the image data to the nearest foreground pixel and at least 
one feature is extracted from the result of the operation. 
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39. A method according to claim 38, wherein at least one 
feature relates to a mean value and/or standard deviation 
and/or the coefficient of variation and/or skewness and/or 
kurtosis of the determined measures for all background 

5 pixels - 

40. A method according to claim 38 or 39, wherein the measure 
for each background pixel is a mean value and/or standard 
deviation and/or the coefficient of variation and/or skewness 
and/or kurtosis and/or maximum distance of the determined 

10 distances in the given directions. 

41. A method according to any of the preceding claims, 
wherein at least one feature related to the bone or to the 
vertebrate, such as age and/ or sex, and/or species, and/or 
race and/or the specific bone considered in the vertebrate, 

15 and/or an estimated Bone Mineral Density of the bone, and/or 
an estimated Bone Mineral Content of the bone, is included in 
the estimation procedure. 

42. A method according to claim 41, wherein the Bone Mineral 
Density is estimated by including data from a reference 

20 object in the exposure of the bone to the electromagnetic 

radiation and where the Bone Mineral Density is estimated by 
correlating the absorption of the electromagnetic radiation 
of the bone and of the reference object. 

43. A method according to any of the preceding claims, 

25 wherein the estimation procedure is based on a statistical 
model, taking into account the correlation structure in the 
data set, so as to assign appropriate weights to the 
significant features in accordance with the predetermined 
relationship . 

30 44 . A method according to any of the preceding claims, 

wherein the reference bone quality is an cUosolute or relative 
bone quality of bones for which image data have been: 
processed in accordance with any of the preceding claims. 
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45. A method according to claim 44, wherein the reference 
bone quality is determined on the basis of measurement of the 
mechanical bone strength, and/or Bone Mineral Density 
measurement, and/ or on the basis of a Bone Mineral Content 

5 measurement, and/or on the basis of a score of the bone 
strength by a skilled radiologist. 

46. A method according to any of the preceding claims, 
wherein the predetermined relationship between the features 
and reference bone quality parameters is estsUalished on the 

10 basis of reference bone quality parameters and features 
extracted according to any of the preceding claims. 

47. A method according to claim 46, wherein the predetermined 
relationship is defined in terms of a model taken from the 
group consisting of: a General Linear Model, a Generalized 

15 Linear Model, an Artificial Neural Network, a Causal 

Probcdailistic Net and Classification And Regression Trees. 

48. A method according to any of the preceding claims, 
wherein the vertebrate is a hxunan. 

49. A method according to any of the preceding claims, 
20 wherein the vertebrate is a horse. 

50. A method according to any of the preceding claims, 
wherein the vertebrate is a large ape, a great ape or an 
anthropoid ape. 

51. A method according to any of the preceding claims, 
25 wherein the vertebrate is a pig. 

52. A method according to any of the preceding claims, 
wherein the vertebrate is a cow. 

53. A method according to any of the preceding claims, 
wherein the bone is taken from the group consisting of 

30 radius, feimir, corpus vertebrae (LI, L2 , L3 , L4 , L5 , Tl, T2 , 
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T3, T4, T5, T6, T7 , T8 . T9 , TIO, Til, T12 , CI, C2 , C3 , C4 , 
C5, C6, C7) , calcaneus, talus, os carpi, metatars, metacarpi, 
falanges, tibia, fibula, patella, ulna, humerus, mandible, 
clavicula, scapula, os coxae, os naviculare, os cuboideum, os 
5 cuneiform I, os cuneiform II, os cuneiform III, os sacrum, os 
coccygis . 

54. A method according to claim 1, wherein procedure (d) is 
used and wherein the features relating to the periodicity of 
the trabecular structure of the part of the bone are derived 
10 from information relating to the height, width and/or 

curvature of peaks of and/or to the total area under peaks of 
the power spectrum of the image data, the smoothed 
periodogram of the image data, or the estimate of the 
parametric estimate of the power spectrum of the image data. 



15 55. A method according to claim 54, wherein the information 
is derived relating to the height, width and/or curvature of 
peaks of the smoothed periodogram of the image data and/or 
total area under the peaks thereof, and wherein the position 
of the peaks are identified in the power spectrum of the 

20 image data. 

56. A method for estimating the bone quality of a vertebrate, 
on the basis of two-dimensional image data comprising 
information relating to the trabecular structure of at least 
a part of a bone of the vertebrate, the image data being data 
25 obtained by exposing at least the part of the bone to 

electromagnetic radiation, the method comprising subjecting 
the image data to a statistical analysis comprising 



a background correction procedure in which low freqpaent 
intensity variations not related to the traibecular 
30 structure of the bone is reduced relative to image data 

related to the traLbecular structure of the part of the 
bone , 
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an image manipulation and feature extraction procedure 
comprising subjecting the image data to at least one of 
the following procedures : 

(a) obtaining an estimate of the parametric estimate 
5 of the power spectrum of the image data and 

extracting features relating to the energy 
distribution of the parametric estimate, 

(b) obtaining, on the basis of image data on which a 
Fourier method has been used to emphasize the 

10 information in the image data relating to the 

trabecular structure, an estimate of a grey- level 
co-occurrence matrix and extracting at least one 
feature on the basis of the estimated co- 
occurrence matrix, 
15 (c) obtaining an estimate of the projected trabecular 

pattern of the image data by using a Fourier 
method to eitphasize the information in the image 
data relating to the trabecular structure and 
subjecting the manipulated image data to a 
20 morphological operation, and extracting features 

relating to the trabecular structure from the 
estimated projected trabecular pattern, 
(d) obtaining, on the basis of a frequency analysis 

the power spectrum of the image data, a smoothed 
25 periodogram of the image data, or a parametric 

estimate of the power spectrtjm of the image data, 
features relating to the periodicity of the 
trabecular structure of the part of the bone, 
and an estimation procedure in which the bone: quality of 
30 the vertebrate is estimated with a Multiple Correlation 

Coefficient better than 0.5 on the basis of the derived 
features and optionally other features related to the 
bone or the vertebrate and a predetermined relationship 
between the features and reference bone quality 
35 parameters. 

57. A method according to claim 56, wherein the bone quality 
of the vertebrate is estimated with a Multiple Correlation 
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Coefficient better than 0.55, such as better than 0.6, 
preferably better than 0.65, such as better than 0.7, 
preferably better than 0.8, such as better than 0.85. 

58. A method for estimating Non-Causal Simultaneous Auto- 
5 Regressive Moving Average models in two or more dimensions," 
the method comprising optimizing a given direct measure of 
the flatness of the residual spectrum of the model. 



59. A method according to claim 58, wherein the optimization 
of the given flatness measure is obtained by 
^° generating a set of initial parameters for the model, 

(b) generating the residual spectrum of the model on the 
basis of the parameters, 

(c) obtaining the measure of the flatness of the residual 
spectrum, 

15 (d) obtaining a new iterate of the parameters on the 

basis of the flatness measure and a search direction 
in parameter space, 
(e) repeating steps (b) - (d) until given stop criterion is 

reached. 

2 0 60. A method according to claims 58 or 5, wherein the measure 
of flatness of the residual spectnam is obtained as the 
entropy of the normalized residual power spectrum, dis- 
regarding the DC value. 
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MORSE estimates vs. true parameters (isotropic case) 
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